%% Estimate IRF Confidence Bounds
% Assume that the data generating process (DGP) is the AR(1) model
% x_t = 1 + 0.9x_{t-2} + u_t,
% where u_t is a series of independent and identically distributed Gaussian variables with mean 0 and variance 1.
% Simulate 500 observations from the model.

rng(1); % For reproducibility
DGP = arima('Constant',1,'AR',{0 0.9},'Variance',1);
y = simulate(DGP,500);

% Explicitly create a diffuse state-space model template for estimation that represents the model.  Fit the model to the data, and return parameter estimates and their corresponding estimated covariance matrix.

A = [0 NaN NaN; 0 1 0; 1 0 0];
B = [NaN; 0; 0];
C = [1 0 0];
D = 0;
Mdl = dssm(A,B,C,D,'StateType',[0 1 0]);
[~,estParams,EstParamCov] = estimate(Mdl,y,abs(randn(3,1)));

% Mdl is an ssm model template for estimation. estParams is a 3-by-1 vector of estimated coefficients. EstParamCov is a 3-by-3 estimated covariance matrix of the coefficient estimates.
% Estimate the IRFs of the state and measurement variables with 95% confidence intervals.

[ResponseY,ResponseX,LowerY,UpperY,LowerX,UpperX] = irf(Mdl,'Params',estParams,...
    'EstParamCov',EstParamCov);

% ResponseY, LowerY, and UpperY are 20-by-1 vectors representing the measurement variable IRF and corresponding lower and upper confidence bounds. ResponseX, LowerX, and UpperX are 20-by-1-by-3 arrays representing the IRF and corresponding lower and upper confidence bounds of the state variables.
% Display a table containing the IRF and confidence bounds of the first state, which represents the AR(2) model.

table(LowerX(:,1,1),ResponseX(:,1,1),UpperX(:,1,1),...
    'VariableNames',["LowerIRFx" "IRFX" "UpperIRFX"])

% The model has only one lag term (lag 2). Therefore, as the shock filters through the system, it impacts the first state variable during odd periods only.
% Copyright 2020 The MathWorks, Inc.