%% Simulate Multiple Response Paths
% Consider this VEC(1) model for three hypothetical response series.
% \begin{array}{rcl}
% \Delta {y_t} &=& c + AB^\prime {y_{t - 1}} + {\Phi _1}\Delta {y_{t - 1}} + {\varepsilon _t}\\
%  & & \\
%  &=& \left[ {\begin{array}{*{20}{c}}
% { - 1}\\
% { - 3}\\
% { - 30}
% \end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
% { - 0.3}&{0.3}\\
% { - 0.2}&{0.1}\\
% { - 1}&0
% \end{array}} \right]\left[ {\begin{array}{*{20}{c}}
% {0.1}&{ - 0.2}&{0.2}\\
% { - 0.7}&{0.5}&{0.2}
% \end{array}} \right]{y_{t - 1}} + \left[ {\begin{array}{*{20}{c}}
% 0&{0.1}&{0.2}\\
% {0.2}&{ - 0.2}&0\\
% {0.7}&{ - 0.2}&{0.3}
% \end{array}} \right]\Delta {y_{t - 1}} + {\varepsilon _t}.
% \end{array}
% The innovations are multivariate Gaussian with a mean of 0 and the covariance matrix
% \Sigma  = \left[
% {\begin{array}{*{20}{c}}{1.3}&{0.4}&{1.6}\\{0.4}&{0.6}&0.7\\{1.6}&0.7&5\end{array}}
% \right].
% Create variables for the parameter values.

Adjustment = [-0.3 0.3; -0.2 0.1; -1 0];
Cointegration = [0.1 -0.7; -0.2 0.5; 0.2 0.2];
ShortRun = {[0. 0.1 0.2; 0.2 -0.2 0; 0.7 -0.2 0.3]};
Constant = [-1; -3; -30];
Trend = [0; 0; 0];
Covariance = [1.3 0.4 1.6; 0.4 0.6 0.7; 1.6 0.7 5];

%
% Create a vecm model object representing the VEC(1) model using the appropriate name-value pair arguments.

Mdl = vecm('Adjustment',Adjustment,'Cointegration',Cointegration,...
    'Constant',Constant,'ShortRun',ShortRun,'Trend',Trend,...
    'Covariance',Covariance);

% Mdl is effectively a fully specified vecm model object. That is, the cointegration constant and linear trend are unknown, but are not needed for simulating observations or forecasting given that the overall constant and trend parameters are known.
% Simulate 1000 paths of 100 observations. Return the innovations (scaled disturbances).
%

numpaths = 1000;
numobs = 100;
rng(1); % For reproducibility
[Y,E] = simulate(Mdl,numobs,'NumPaths',numpaths);

% Y is a 100-by-3-by-1000 matrix of simulated responses. E is a matrix whose dimensions correspond to the dimensions of Y, but represents the simulated, scaled disturbances. Columns correspond to the response variable names Mdl.SeriesNames.
%
% For each time point, compute the mean vector of the simulated responses among all paths.

MeanSim = mean(Y,3);

% MeanSim is a 100-by-7 matrix containing the average of the simulated responses at each time point.
%
% Plot the simulated responses and their averages.

figure;
for j = 1:Mdl.NumSeries
    subplot(2,2,j)
    plot(squeeze(Y(:,j,:)),'Color',[0.8,0.8,0.8])
    title(Mdl.SeriesNames{j});
    hold on
    plot(MeanSim(:,j));
    xlabel('Time index')
    hold off
end

% Copyright 2012 The MathWorks, Inc.