%% Forecasting
% A TSP and a DSP produce different forecasts. Basically, shocks to a TSP return to the trend line c+\delta t as time increases. In contrast, shocks to a DSP might be persistent over time.
% For example, consider the simple trend-stationary model
% y_{1,t}=0.9y_{1,t-1}+0.02t+\varepsilon_{1,t}
% and the difference-stationary model
% y_{2,t}=0.2+y_{2,t-1}+\varepsilon_{2,t}.
% In these models, \varepsilon_{1,t} and \varepsilon_{2,t} are independent innovation processes. For this example, the innovations are independent and distributed N(0,1).
% Both processes grow at rate 0.2. To calculate the growth rate for the TSP, which has a linear term 0\ldotp 02t, set \varepsilon_{1,t} =0. Then, solve the model y_{1,t} =c+\delta t for c and \delta.
% c+\delta t=0.9(c+\delta(t-1))+0.02t.
% The solution is c=-1\ldotp 8, \delta =0\ldotp 2.
% A plot for t = 1:1000 shows the TSP stays very close to the trend line, while the DSP has persistent deviations away from the trend line.

T = 1000;   % Sample size
t = (1:T)'; % Period vector
rng(5);     % For reproducibility

randm = randn(T,2); % Innovations
y = zeros(T,2);     % Columns of y are data series

% Build trend stationary series
y(:,1) = .02*t + randm(:,1); 
for ii = 2:T
    y(ii,1) = y(ii,1) + y(ii-1,1)*.9;
end

% Build difference stationary series
y(:,2) = .2 + randm(:,2); 
y(:,2) = cumsum(y(:,2)); 

figure
plot(y(:,1),'b')
hold on
plot(y(:,2),'g')
plot((1:T)*0.2,'k--')
legend('Trend Stationary','Difference Stationary',...
    'Trend Line','Location','NorthWest')
hold off

%
% Forecasts based on the two series are different. To see this difference, plot the predicted behavior of the two series using docid:econ_ug#bvn5drn, docid:econ_ug#bvn5d8w, and docid:econ_ug#bvn64fj. The following plot shows the last 100 data points in the two series and predictions of the next 100 points, including confidence bounds.

AR = {[NaN 0; 0 NaN]}; % Independent response series
trend = [NaN; 0];      % Linear trend in first series only
Mdl = varm('AR',AR,'Trend',trend);

EstMdl = estimate(Mdl,y);
EstMdl.SeriesNames = ["Trend stationary" "Difference stationary"];

[ynew,ycov] = forecast(EstMdl,100,y);
% This generates predictions for 100 time steps

seY = sqrt(diag(EstMdl.Covariance))'; % Extract standard deviations of y
CIY = zeros([size(y) 2]);             % In-sample intervals
CIY(:,:,1) = y - seY;
CIY(:,:,2) = y + seY;

extractFSE = cellfun(@(x)sqrt(diag(x))',ycov,'UniformOutput',false);
seYNew = cell2mat(extractFSE);
CIYNew = zeros([size(ynew) 2]); % Forecast intervals
CIYNew(:,:,1) = ynew - seYNew;
CIYNew(:,:,2) = ynew + seYNew;

tx = (T-100:T+100);
hs = 1:2;
figure;
for j = 1:Mdl.NumSeries
    hs(j) = subplot(2,1,j);
    hold on;
    h1 = plot(tx,tx*0.2,'k--');
    axis tight;
    ha = gca;
    h2 = plot(tx,[y(end-100:end,j); ynew(:,j)]);
    h3 = plot(tx(1:101),squeeze(CIY(end-100:end,j,:)),'r:');
    plot(tx(102:end),squeeze(CIYNew(:,j,:)),'r:');
    h4 = fill([tx(102) ha.XLim([2 2]) tx(102)],ha.YLim([1 1 2 2]),[0.7 0.7 0.7],...
        'FaceAlpha',0.1,'EdgeColor','none');
    title(EstMdl.SeriesNames{j});
    hold off;
end
legend(hs(1),[h1 h2 h3(1) h4],...
    {'Trend','Process','Interval estimate','Forecast horizon'},'Location','Best');

% Examine the fitted parameters by passing the estimated model to summarize, and you find estimate did an excellent job.
% The TSP has confidence intervals that do not grow with time, whereas the DSP has confidence intervals that grow. Furthermore, the TSP goes to the trend line quickly, while the DSP does not tend towards the trend line y=0\ldotp 2t asymptotically.
% Copyright 2020 The MathWorks, Inc.