%% Simulate AR Process
% This example shows how to simulate sample paths from a stationary AR(2) process without specifying presample observations.
%% Create Model
% Specify the AR(2) model
% {y_t} = 0.5 + 0.7{y_{t - 1}} + 0.25{y_{t - 2}} + {\varepsilon _t},
% where the innovation process is Gaussian with variance 0.1.

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

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

rng(5)
Y = simulate(Mdl,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
% \mu = \frac{c}{{\left( {1 - {\phi _1} - {\phi _2}} \right)}}.
% The unconditional variance of the process is
% \sigma^2 = \frac{{(1 - {\phi _2})}}{{(1 + {\phi _2})}}\frac{{\sigma _\varepsilon ^2}}{{{{(1 - {\phi _2})}^2} - \phi _1^2}}.
% Compute the unconditional mean and variance of the process.

c = Mdl.Constant;
phi = Mdl.AR;
sigmaEps2 = Mdl.Variance;
mu = c/(1-phi{1}-phi{2})
sigma2 = (1-phi{2})/(1+phi{2})*sigmaEps2/((1-phi{2})^2 - phi{1}^2)

%
%% Generate Many Sample Paths
% Generate 1000 sample paths, each with 50 observations.

Y = simulate(Mdl,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,sigma2*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.
% Around observation 50, the simulated variance approaches the theoretical variance.
%
%% Oversample to Reduce Transient Effects
% 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.

Y = simulate(Mdl,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,sigma2*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.
% Copyright 2015 The MathWorks, Inc.