Documentation

## Simulate Stationary Processes

### Simulate an AR Process

This example shows how to simulate sample paths from a stationary AR(2) process without specifying presample observations.

Step 1. Specify a model.

Specify the AR(2) model

`${y}_{t}=0.5+0.7{y}_{t-1}+0.25{y}_{t-2}+{\epsilon }_{t},$`

where the innovation process is Gaussian with variance 0.1.

`model = arima('Constant',0.5,'AR',{0.7,0.25},'Variance',.1);`

Step 2. Generate one sample path.

Generate one sample path (with 50 observations) from the specified model, and plot.

```rng('default') Y = simulate(model,50); figure plot(Y) xlim([0,50]) title('Simulated AR(2) Process')``` Because presample data was not specified, `simulate` sets the two required presample observations equal to the unconditional mean of the process,

`$\frac{c}{\left(1-{\varphi }_{1}-{\varphi }_{2}\right)}=\frac{0.5}{\left(1-0.7-0.25\right)}=10.$`

Step 3. Generate many sample paths.

Generate 1000 sample paths, each with 50 observations.

```rng('default') Y = simulate(model,50,'NumPaths',1000); figure subplot(2,1,1) plot(Y,'Color',[.85,.85,.85]) title('Simulated AR(2) Process') hold on h=plot(mean(Y,2),'k','LineWidth',2); legend(h,'Simulation Mean','Location','NorthWest') hold off subplot(2,1,2) plot(var(Y,0,2),'r','LineWidth',2) title('Process Variance') hold on plot(1:50,.83*ones(50,1),'k--','LineWidth',1.5) legend('Simulation','Theoretical',... 'Location','SouthEast') hold off``` The simulation mean is constant over time. This is consistent with the definition of a stationary process. The process variance is not constant over time, however. There are transient effects at the beginning of the simulation due to the absence of presample data.

The simulated variance approaches the theoretical variance,

`$\frac{\left(1-{\varphi }_{2}\right)}{\left(1+{\varphi }_{2}\right)}\frac{{\sigma }_{\epsilon }^{2}}{{\left(1-{\varphi }_{2}\right)}^{2}-{\varphi }_{1}^{2}}=\frac{\left(1-.25\right)}{\left(1+.25\right)}\frac{0.1}{{\left(1-.25\right)}^{2}-{0.7}^{2}}=0.83,$`

by around the 50th observation.

Step 4. Oversample the process.

To reduce transient effects, one option is to oversample the process. For example, to sample 50 observations, you can generate paths with more than 50 observations, and discard all but the last 50 observations as burn-in. Here, simulate paths of length 150, and discard the first 100 observations.

```rng('default') Y = simulate(model,150,'NumPaths',1000); Y = Y(101:end,:); figure subplot(2,1,1) plot(Y,'Color',[.85,.85,.85]) title('Simulated AR(2) Process') hold on h=plot(mean(Y,2),'k','LineWidth',2); legend(h,'Simulation Mean','Location','NorthWest') hold off subplot(2,1,2) plot(var(Y,0,2),'r','LineWidth',2) xlim([0,50]) title('Process Variance') hold on plot(1:50,.83*ones(50,1),'k--','LineWidth',1.5) legend('Simulation','Theoretical',... 'Location','SouthEast') hold off``` The realizations now look like draws from a stationary stochastic process. The simulation variance fluctuates (due to Monte Carlo error) around the theoretical variance.

### Simulate an MA Process

This example shows how to simulate sample paths from a stationary MA(12) process without specifying presample observations.

Step 1. Specify a model.

Specify the MA(12) model

`${y}_{t}=0.5+{\epsilon }_{t}+0.8{\epsilon }_{t-1}+0.2{\epsilon }_{t-12},$`

where the innovation distribution is Gaussian with variance 0.2.

```model = arima('Constant',0.5,'MA',{0.8,0.2},... 'MALags',[1,12],'Variance',0.2);```

Step 2. Generate sample paths.

Generate 200 sample paths, each with 60 observations.

```rng('default') Y = simulate(model,60,'NumPaths',200); figure plot(Y,'Color',[.85,.85,.85]) hold on h = plot(mean(Y,2),'k','LineWidth',2)```
```h = Line with properties: Color: [0 0 0] LineStyle: '-' LineWidth: 2 Marker: 'none' MarkerSize: 6 MarkerFaceColor: 'none' XData: [1x60 double] YData: [1x60 double] ZData: [1x0 double] Show all properties ```
```legend(h,'Simulation Mean','Location','NorthWest') title('MA(12) Process') hold off``` For an MA process, the constant term is the unconditional mean. The simulation mean is around 0.5, as expected.

Step 3. Plot the simulation variance.

The unconditional variance for the model is

`$\left(1+{\theta }_{1}^{2}+{\theta }_{12}^{2}\right){\sigma }_{\epsilon }^{2}=\left(1+0.{8}^{2}+0.{2}^{2}\right)×0.2=0.336.$`

Because the model is stationary, the unconditional variance should be constant across all times. Plot the simulation variance, and compare it to the theoretical variance.

```figure plot(var(Y,0,2),'Color',[.75,.75,.75],'LineWidth',1.5) xlim([0,60]) title('Unconditional Variance') hold on plot(1:60,.336*ones(60,1),'k--','LineWidth',2) legend('Simulation','Theoretical',... 'Location','SouthEast') hold off``` There appears to be a short burn-in period at the beginning of the simulation. During this time, the simulation variance is lower than expected. Afterwards, the simulation variance fluctuates around the theoretical variance.

Step 4. Generate more sample paths.

Simulate 10,000 paths from the model, each with length 1000. Look at the simulation variance.

```rng('default') YM = simulate(model,1000,'NumPaths',10000); figure plot(var(YM,0,2),'Color',[.75,.75,.75],'LineWidth',1.5) ylim([0.3,0.36]) title('Unconditional Variance') hold on plot(1:1000,.336*ones(1000,1),'k--','LineWidth',2) legend('Simulation','Theoretical',... 'Location','SouthEast') hold off``` The Monte Carlo error is reduced when more realizations are generated. There is much less variability in the simulation variance, which tightly fluctuates around the theoretical variance.