%% IRF of Seasonal AR Model
% This example shows how to plot and return the IRF of a seasonal AR model by using impulse and armairf. Also, the example shows how to prepare LagOp lag operator polynomials as inputs to armairf.
% The equation of an SAR\left(p,0,0\right)\times \left(p_s ,0,0{\left.\right)}_s \right. model is
% y_t=c+\Phi(L)^{-1}\phi(L)^{-1}\varepsilon_t,
% where:
% \phi \left(L\right) is the p-degree AR lag operator polynomial 1-\phi_1 L-\ldotp \ldotp \ldotp -\phi_p L^p.
% \Phi \left(L\right) is the p_s-degree seasonal AR lag operator polynomial 1-\Phi_{p_1 } L^{p_1 } -\ldotp \ldotp \ldotp -\Phi_{p_s } L^{p_s }.
% Like a pure AR process, an SAR process is stationary when the product \phi \left(L\right)\Phi \left(L\right) is stable. In this case, the infinite-degree inverse polynomial \psi (L) = \Phi(L)^{-1}\phi {(L)^{ - 1}} has absolutely summable coefficients, and the IRF decays to zero.
%
%% Plot IRF Using impulse
% Create a quarterly SAR\left(1,0,0\right)\times {\left(4,0,0\right)}_4 model with coefficients \phi_1 =0\ldotp 5 and \Phi_4 =-0\ldotp 4, a model constant of 0, and an innovation variance of 1.

ar = 0.5;
sar = -0.4;
sarlags = 4;
Mdl = arima('AR',ar,'SAR',sar,'SARLags',sarlags,...
    'Constant',0,'Variance',1)

% Plot the IRF of the SAR model for 17 quarters, from quarters 0 through 16.

numObs = 17;
impulse(Mdl,numObs)

% The IRF decays in a sinusoidal pattern.
% Return the IRF.

irfIMPULSE = impulse(Mdl,numObs);

%
%% Plot IRF Using armirf
% armairf accepts one overall AR polynomial. Therefore, you must multiply all AR and differencing lag operator polynomials present in the model before calling armairf.
% Create lag operator polynomials for the AR and SAR polynomials. For each polynomial:
% Include the lag 0 term, which has a coefficient of 1.
% Negate the coefficients to express the polynomials in lag operator notation, with all AR polynomials on the left side of the equation.

ARLOP = LagOp([1 -ar],'Lags',[0 1])
MALOP = LagOp([1 -sar],'Lags',[0 sarlags])

% Multiply the polynomials.

ARProdLOP = ARLOP*MALOP

% ARProdLOP is a LagOp object representing the product of the AR and SAR polynomials of the SAR model.
% Plot and return the IRF by passing ARProdLOP as the AR polynomial (first input). Specify an empty array for the MA polynomial (second input). To plot the IRF, also return the plot handle.

[irfARMAIRF,h] = armairf(ARProdLOP,[],'NumObs',numObs);

% Compare the IRFs.

periods = (0:(numObs - 1))';
table(periods,irfIMPULSE,irfARMAIRF)

% The IRFs returned by the two functions appear equivalent.
% Copyright 2021 The MathWorks, Inc.