%% Simulate MA Process
% This example shows how to simulate sample paths from a stationary MA(12) process without specifying presample observations.
%% Create Model
% Specify the MA(12) model
% {y_t} = 0.5 + {\varepsilon _t} + 0.8{\varepsilon _{t - 1}} + 0.2{\varepsilon _{t - 12}},
% where the innovation distribution is Gaussian with variance 0.2.

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

%
%% Generate Sample Paths
% Generate 200 sample paths, each with 60 observations.

rng(5)
Y = simulate(Mdl,60,'NumPaths',200);

figure
plot(Y,'Color',[.85,.85,.85])
hold on
h = plot(mean(Y,2),'k','LineWidth',2);
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 approximately c = 0.5.
%
%% Plot Simulation Variance
% The unconditional variance of the process is
% \sigma^2 = (1 + \theta _1^2 + \theta _{12}^2)\sigma _\varepsilon ^2.
% Compute the unconditional variance.

theta = cell2mat(Mdl.MA);
sigmaEps2 = Mdl.Variance;
sigma2 = (1+sum(theta.^2))*sigmaEps2

% 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,sigma2*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.
%
%% Generate Many Sample Paths
% Simulate 10,000 paths from the model, each with length 1000. Plot the simulation variance.

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