%% Forecast Model with VARX Submodels
% Compute optimal and estimated forecasts and corresponding forecast error covariance matrices from a three-state Markov-switching dynamic regression model for a 2-D VARX response process. This example uses arbitrary parameter values for the DGP.
%
%% Create Fully Specified Model for DGP
% Create a three-state discrete-time Markov chain model for the switching mechanism.

P = [10 1 1; 1 10 1; 1 1 10];
mc = dtmc(P);

% mc is a fully specified dtmc object.  dtmc normalizes the rows of P so that they sum to 1.
% For each regime, use varm to create a VAR model that describes the response process within the regime. Specify all parameter values.

% Constants
C1 = [1;-1];
C2 = [2;-2];
C3 = [3;-3];

% Autoregression coefficients
AR1 = {};                            
AR2 = {[0.5 0.1; 0.5 0.5]};          
AR3 = {[0.25 0; 0 0] [0 0; 0.25 0]}; 

% Regression coefficients
Beta1 = [1;-1];           
Beta2 = [2 2;-2 -2];      
Beta3 = [3 3 3;-3 -3 -3]; 

% Innovations covariances
Sigma1 = [1 -0.1; -0.1 1];
Sigma2 = [2 -0.2; -0.2 2];
Sigma3 = [3 -0.3; -0.3 3];

% VARX submodels
mdl1 = varm('Constant',C1,'AR',AR1,'Beta',Beta1,'Covariance',Sigma1);
mdl2 = varm('Constant',C2,'AR',AR2,'Beta',Beta2,'Covariance',Sigma2);
mdl3 = varm('Constant',C3,'AR',AR3,'Beta',Beta3,'Covariance',Sigma3);
mdl = [mdl1; mdl2; mdl3];

% mdl contains three fully specified varm model objects.
% For the DGP, create a fully specified Markov-switching dynamic regression model from the switching mechanism mc and the submodels mdl.

Mdl = msVAR(mc,mdl);

% Mdl is a fully specified msVAR model.
%
%% Forecast Model Ignoring Regression Component
% If you do not supply exogenous data, simulate and forecast ignore the regression components in the submodels. forecast requires enough data before the forecast horizon to initialize the model. Simulate 120 observations from the DGP.

rng('default'); % For reproducibility
Y = simulate(Mdl,120);

% Y is a 120-by-2 matrix of simulated responses.  Rows correspond to time points, and columns correspond to variables in the system.
% Treat the first 100 observations of the simulated response data as the presample for the forecast, and treat the last 20 observations as a holdout sample.

idx0 = 1:100;
idx1 = 101:120;
Y0 = Y(idx0,:); % Forecast sample
Y1 = Y(idx1,:); % Holdout sample

% Compute 1- through 20-step-ahead optimal and estimated point forecasts from the model. Compute forecast error covariance matrices corresponding to the estimated forecasts.

YF1 = forecast(Mdl,Y0,20); 
[YF2,EstCov] = forecast(Mdl,Y0,20);

% YF1 and YF2 are 20-by-2 matrices of optimal and estimated forecasts, respectively. EstCov is a 2-by-2-by-20 array of forecast error covariances.
% Extract the forecast error variances of each response for each period in the forecast horizon.

estVar1(:) = EstCov(1,1,:);
estVar2(:) = EstCov(2,2,:);

% estVar1 and estVar2 are 1-by-20 vectors of forecast error variances.
% Plot the data, forecasts, and 95% forecast intervals of each variable on separate subplots.

figure

subplot(2,1,1)
hold on
plot(idx0,Y0(:,1),'b');
h = plot(idx1,Y1(:,1),'b--');
h1 = plot(idx1,YF1(:,1),'r');
h2 = plot(idx1,YF2(:,1),'m');
ciu = YF2(:,1) + 1.96*sqrt(estVar1');   % Upper 95% confidence level
cil = YF2(:,1) - 1.96*sqrt(estVar1');   % Lower 95% confidence level
plot(idx1,ciu,'m-.'); 
plot(idx1,cil,'m-.'); 
yfill = [ylim,fliplr(ylim)];
xfill = [idx0(end) idx0(end) idx1(end) idx1(end)];
fill(xfill,yfill,'k','FaceAlpha',0.05)
legend([h h1 h2],["Actual" "Optimal" "Estimated"],...
    'Location',"NorthWest")
title("Point and Interval Forecasts: Series 1")
hold off

subplot(2,1,2)
hold on
plot(idx0,Y0(:,2),'b');
h = plot(idx1,Y1(:,2),'b--');
h1 = plot(idx1,YF1(:,2),'r');
h2 = plot(idx1,YF2(:,2),'m');
ciu = YF2(:,2) + 1.96*sqrt(estVar2');   % Upper 95% confidence level
cil = YF2(:,2) - 1.96*sqrt(estVar2');   % Lower 95% confidence level
plot(idx1,ciu,'m-.'); 
plot(idx1,cil,'m-.'); 
yfill = [ylim,fliplr(ylim)];
xfill = [idx0(end) idx0(end) idx1(end) idx1(end)];
fill(xfill,yfill,'k','FaceAlpha',0.05)
legend([h h1 h2],["Actual" "Optimal" "Estimated"],...
    'Location',"NorthWest")
title("Point and Interval Forecasts: Series 2")
hold off

%
%% Forecast Model Including Regression Component
% Simulate exogenous data for the three regressors by generating 120 random observations from the 3-D standard Gaussian distribution.

X = randn(120,3);

% Simulate 120 observations from the DGP. Specify the exogenous data for the regression component.

rng('default') 
Y = simulate(Mdl,120,'X',X);

% Treat the first 100 observations of the simulated response and exogenous data as the presample for the forecast, and treat the last 20 observations as a holdout sample.

idx0 = 1:100;
idx1 = 101:120;
Y0 = Y(idx0,:);
Y1 = Y(idx1,:);
X1 = X(idx1,:);

% Compute 1- through 20-step-ahead optimal and estimated point forecasts from the model. Compute forecast error covariance matrices corresponding to the estimated forecasts. Specify the forecast-period exogenous data for the regression component.

YF1 = forecast(Mdl,Y0,20,'X',X1); 
[YF2,EstCov] = forecast(Mdl,Y0,20,'X',X1);
estVar1(:) = EstCov(1,1,:);
estVar2(:) = EstCov(2,2,:);

% Plot the data, forecasts, and 95% forecast intervals of each variable on separate subplots.

figure

subplot(2,1,1)
hold on
plot(idx0,Y0(:,1),'b');
h = plot(idx1,Y1(:,1),'b--');
h1 = plot(idx1,YF1(:,1),'r');
h2 = plot(idx1,YF2(:,1),'m');
ciu = YF2(:,1) + 1.96*sqrt(estVar1');   % Upper 95% confidence level
cil = YF2(:,1) - 1.96*sqrt(estVar1');   % Lower 95% confidence level
plot(idx1,ciu,'m-.');
plot(idx1,cil,'m-.');
yfill = [ylim,fliplr(ylim)];
xfill = [idx0(end) idx0(end) idx1(end) idx1(end)];
fill(xfill,yfill,'k','FaceAlpha',0.05)
legend([h h1 h2],["Actual" "Optimal" "Estimated"],...
    'Location',"NorthWest")
title("Point and Interval Forecasts: Series 1")
hold off

subplot(2,1,2)
hold on
plot(idx0,Y0(:,2),'b');
h = plot(idx1,Y1(:,2),'b--');
h1 = plot(idx1,YF1(:,2),'r');
h2 = plot(idx1,YF2(:,2),'m');
ciu = YF2(:,2) + 1.96*sqrt(estVar2');   % Upper 95% confidence level
cil = YF2(:,2) - 1.96*sqrt(estVar2');   % Lower 95% confidence level
plot(idx1,ciu,'m-.'); 
plot(idx1,cil,'m-.'); 
yfill = [ylim,fliplr(ylim)];
xfill = [idx0(end) idx0(end) idx1(end) idx1(end)];
fill(xfill,yfill,'k','FaceAlpha',0.05)
legend([h h1 h2],["Actual" "Optimal" "Estimated"],...
    'Location',"NorthWest")
title("Point and Interval Forecasts: Series 2")
hold off

% Copyright 2019 The MathWorks, Inc.
%