Main Content

irf

Impulse response function (IRF) of state-space model

Since R2020b

Description

irf returns a numeric array representing the IRFs of the state and measurement variables in a state-space model. To plot the IRFs instead, use irfplot. Other state-space model tools to characterize the dynamics of a specified system include:

  • The forecast error variance decomposition (FEVD), computed by fevd, provides information about the relative importance of each state disturbance in affecting the forecast error variance of all measurement variables in the system.

  • Model-implied temporal correlations, computed by corr for a standard state-space model, measure the association between present and past state or measurement variables, as prescribed by the form of the model.

Fully Specified State-Space Model

example

ResponseY = irf(Mdl) returns the IRF, or dynamic response, of each measurement variable ResponseY of the fully specified state-space model Mdl, such as an estimated model.

example

ResponseY = irf(Mdl,Name,Value) uses additional options specified by one or more name-value pair arguments. For example, 'NumPeriods',10,'Cumulative',true specifies a 10-period cumulative IRF starting at time 1, during which irf applies the shock to a state-disturbance variable in the system, and ending at period 10.

example

[ResponseY,ResponseX] = irf(___) also returns the IRF of each state variable ResponseX, using any of the input argument combinations in the previous syntaxes.

Partially Specified State-Space Model and Confidence Interval Estimation

example

[ResponseY,ResponseX] = irf(___,'Params',estParams) returns the IRFs of all variables of the partially specified state-space model Mdl. estParams specifies estimates of all unknown parameters in the model.

example

[ResponseY,ResponseX,LowerY,UpperY,LowerX,UpperX] = irf(___,'Params',estParams,'EstParamCov',EstParamCov) also returns, for each period, the lower and upper 95% Monte Carlo confidence bounds of each measurement variable IRF ([LowerY,UpperY]) and each state variable IRF ([LowerX,UpperX]). EstParamCov specifies the estimated covariance matrix of the parameter estimates, as returned by the estimate function, and is required for confidence interval estimation.

Examples

collapse all

Explicitly create the state-space model

xt=0.5xt-1+0.2utyt=2xt+0.01εt.

A = 0.5;
B = 0.2;
C = 2;
D = 0.01;
Mdl = ssm(A,B,C,D)
Mdl = 
State-space model type: ssm

State vector length: 1
Observation vector length: 1
State disturbance vector length: 1
Observation innovation vector length: 1
Sample size supported by model: Unlimited

State variables: x1, x2,...
State disturbances: u1, u2,...
Observation series: y1, y2,...
Observation innovations: e1, e2,...

State equation:
x1(t) = (0.50)x1(t-1) + (0.20)u1(t)

Observation equation:
y1(t) = (2)x1(t) + (0.01)e1(t)

Initial state distribution:

Initial state means
 x1 
  0 

Initial state covariance matrix
     x1   
 x1  0.05 

State types
     x1     
 Stationary 

Mdl is a ssm model object. Because all parameters have known values, the object is fully specified.

Compute the IRF of the measurement variable.

responseY = irf(Mdl)
responseY = 20×1

    0.4000
    0.2000
    0.1000
    0.0500
    0.0250
    0.0125
    0.0063
    0.0031
    0.0016
    0.0008
      ⋮

responseY is a 20-by-1 vector representing the 20-period IRF of the measurement variable yt. responseY(5) is 0.0250, which means that the response of yt at period 5, to a unit shock to the state disturbance ut at period 1, is 0.0250.

Explicitly create the multivariate diffuse state-space model

x1,t=x1,t-1+0.2u1,tx2,t=x1,t-1+0.3x2,t-1+u2,ty1,t=x1,t+ε1,ty2,t=x1,t+x2,t+ε2,t.

A = [1 0; 1 0.3];
B = [0.2 0; 0 1];
C = [1 0; 1 1];
D = eye(2);
Mdl = dssm(A,B,C,D,'StateType',[2 2])
Mdl = 
State-space model type: dssm

State vector length: 2
Observation vector length: 2
State disturbance vector length: 2
Observation innovation vector length: 2
Sample size supported by model: Unlimited

State variables: x1, x2,...
State disturbances: u1, u2,...
Observation series: y1, y2,...
Observation innovations: e1, e2,...

State equations:
x1(t) = x1(t-1) + (0.20)u1(t)
x2(t) = x1(t-1) + (0.30)x2(t-1) + u2(t)

Observation equations:
y1(t) = x1(t) + e1(t)
y2(t) = x1(t) + x2(t) + e2(t)

Initial state distribution:

Initial state means
 x1  x2 
  0   0 

Initial state covariance matrix
     x1   x2  
 x1  Inf  0   
 x2  0    Inf 

State types
    x1       x2   
 Diffuse  Diffuse 

Mdl is a dssm model object.

Compute the 10-period IRFs of the measurement variables.

ResponseY = irf(Mdl,'NumPeriods',10);

ResponseY is a 10-by-2-by-2 array representing the 10-period IRFs of the measurement variables. For example, ResponseY(:,1,2) is the IRF of y2,t as a result of a shock applied to u1,t.

ResponseY(:,1,2)
ans = 10×1

    0.2000
    0.4000
    0.4600
    0.4780
    0.4834
    0.4850
    0.4855
    0.4857
    0.4857
    0.4857

Simulate data from a known model, fit the data to a state-space model, and then estimate the cumulative IRFs of the state variables.

Assume that the data generating process (DGP) is the AR(1) model

xt=1+0.9xt-2+ut,

where ut 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 state-space model template for estimation that represents the model

xt=c+ϕxt-2+ηutyt=xt.

A = [0 NaN NaN; 0 1 0; 1 0 0];
B = [NaN; 0; 0];
C = [1 0 0];
D = 0;
Mdl = ssm(A,B,C,D,'StateType',[0 1 0]);

Fit the model template to the data. Specify a set of positive, random standard Gaussian starting values for the three model parameters.

EstMdl = estimate(Mdl,y,abs(randn(3,1)));
Method: Maximum likelihood (fminunc)
Sample size: 500
Logarithmic  likelihood:     -2085.74
Akaike   info criterion:      4177.49
Bayesian info criterion:      4190.13
      |     Coeff       Std Err   t Stat     Prob  
---------------------------------------------------
 c(1) |  0.36553       0.07967    4.58829  0.00000 
 c(2) |  0.70179       0.00738   95.13852   0      
 c(3) |  1.16649       0.02236   52.16929   0      
      |                                            
      |   Final State   Std Dev    t Stat    Prob  
 x(1) | 10.72536        0          Inf      0      
 x(2) |   1             0          Inf      0      
 x(3) |  6.66084        0          Inf      0      

EstMdl is a fully specified ssm model object.

Estimate the cumulative IRFs of the state and measurement variables.

[ResponseY,ResponseX] = irf(EstMdl,'Cumulative',true);

ResponseY is a 20-by-1 vector representing the measurement variable IRF. ResponseX is a 20-by-1-by-3 array representing the IRF of the state variables.

Display the IRF of xt, which is the first state variable in the system x1,t.

irfx = ResponseX(:,:,1)
irfx = 20×1

    1.1665
    1.1665
    1.9851
    1.9851
    2.5596
    2.5596
    2.9628
    2.9628
    3.2458
    3.2458
      ⋮

Verify that, because yt=xt, ResponseY = ResponseX(:,:,1).

ver1 = sum(abs(ResponseY - ResponseX(:,:,1)))
ver1 = 0

Verify that, because x1,t-1=x3,t, ResponseX(1:(end-2),1,1) = ResponseX(2:(end-1),:,3).

ver2 = sum(abs(ResponseX(1:(end-2),:,1) - ResponseX(2:(end-1),:,3)))
ver2 = 0

Simulate data from a time-varying state-space model, fit a model to the data, and then estimate the time-varying IRF.

Consider the DGP represented by the system

xt={0.75xt-1+ut;t<11-0.1xt-1+3ut;t11yt=1.5xt+2εt.

Write a function that specifies how the parameters params map to the state-space model matrices, the initial state moments, and the state types. Save this code as a file named timeVariantAR1ParamMap.m on your MATLAB® path. Alternatively, open the example to access the function.

type timeVariantAR1ParamMap.m
% Copyright 2020 The MathWorks, Inc.

function [A,B,C,D] = timeVariantAR1ParamMap(params)
% Time-varying state-space model parameter mapping function example. This
% function maps the vector params to the state-space matrices (A, B, C, and
% D). From periods 1 through 10, the state model is an AR(1)model, and from
% periods 11 through 20, the state model is possibly a different AR(1)
% model. The measurement equation is the same throughout the time span.
    A1 = {params(1)};
    A2 = {params(2)};
    varu1 = exp(params(3));  % Positive variance constraints
    varu2 = exp(params(4));
    B1 = {sqrt(varu1)}; 
    B2 = {sqrt(varu2)};
    C = params(5);
    vare1 = exp(params(6));
    D = sqrt(vare1);
    A = [repmat(A1,10,1); repmat(A2,10,1)];
    B = [repmat(B1,10,1); repmat(B2,10,1)];
end

Implicitly create a partially specified state-space model representing the DGP. For this example, fix the measurement-sensitivity coefficient C to 1.5.

C = 1.5;
fixCParamMap = @(x)timeVariantAR1ParamMap([x(1:4), C, x(5)]);
DGP = ssm(fixCParamMap);

Simulate 20 observations from the DGP. Because DGP is partially specified, pass the true parameter values to simulate by using the 'Params' name-value pair argument.

rng(10) % For reproducibility
A1 = 0.75;
A2 = -0.1; 
B1 = 1;
B2 = 3;
D = 2;
trueParams = [A1 A2 2*log(B1) 2*log(B2) 2*log(D)]; % Transform variances for parameter map
y = simulate(DGP,20,'Params',trueParams);

y is a 20-by-1 vector of simulated measurements yt from the DGP.

Because DGP is a partially specified, implicit model object, its parameters are unknown. Therefore, it can serve as a model template for estimation.

Fit the model to the simulated data. Specify random standard Gaussian draws for the initial parameter values. Return the parameter estimates.

[~,estParams] = estimate(DGP,y,randn(1,5),'Display','off')
estParams = 1×5

    0.6164   -0.1665    0.0135    1.6803   -1.5855

estParams is a 1-by-5 vector of parameter estimates. The output argument list of the parameter mapping function determines the order of the estimates: A{1}, A{2}, B{1}, B{2}, and D.

Estimate the IRFs of the measurement and state variables by supplying DGP (not the estimated model) and the estimated parameters using the 'Params' name-value pair argument.

[responseY,responseX] = irf(DGP,'Params',estParams);
table(responseY,responseX)
ans=20×2 table
     responseY      responseX 
    ___________    ___________

         1.5101         1.0068
        0.93091         0.6206
        0.57385        0.38257
        0.35374        0.23583
        0.21806        0.14537
        0.13442       0.089615
       0.082863       0.055242
        0.05108       0.034054
       0.031488       0.020992
       0.019411        0.01294
     -0.0032311     -0.0021541
     0.00053785     0.00035857
    -8.9531e-05    -5.9687e-05
     1.4903e-05     9.9356e-06
    -2.4808e-06    -1.6539e-06
     4.1296e-07     2.7531e-07
      ⋮

responseY and responseX are time-varying IRFs. The first 10 periods correspond to the IRF of the first state equation. During period 11, the remainder of the shock transfers to the second state equation and filters through that system until it diminishes.

Assume that the data generating process (DGP) is the AR(1) model

xt=1+0.9xt-2+ut,

where ut 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)));
Method: Maximum likelihood (fminunc)
Effective Sample size:            500
Logarithmic  likelihood:     -2085.74
Akaike   info criterion:      4177.49
Bayesian info criterion:      4190.13
      |     Coeff       Std Err   t Stat     Prob  
---------------------------------------------------
 c(1) |  0.36553       0.07967    4.58829  0.00000 
 c(2) |  0.70179       0.00738   95.13852   0      
 c(3) |  1.16649       0.02236   52.16929   0      
      |                                            
      |   Final State   Std Dev    t Stat    Prob  
 x(1) | 10.72536        0          Inf      0      
 x(2) |   1             0          Inf      0      
 x(3) |  6.66084        0          Inf      0      

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"])
ans=20×3 table
    LowerIRFx      IRFX      UpperIRFX
    _________    ________    _________

      1.1214       1.1665       1.209 
           0            0           0 
     0.78826      0.81864     0.84833 
           0            0           0 
     0.54845      0.57452     0.60214 
           0            0           0 
     0.37964      0.40319     0.42929 
           0            0           0 
      0.2609      0.28296     0.30597 
           0            0           0 
     0.17908      0.19858     0.21954 
           0            0           0 
     0.12339      0.13936     0.15655 
           0            0           0 
    0.084751     0.097803     0.11184 
           0            0           0 
      ⋮

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.

Input Arguments

collapse all

State-space model, specified as an ssm model object returned by ssm or its estimate function, or a dssm model object returned by dssm or its estimate function.

If Mdl is partially specified (that is, it contains unknown parameters), specify estimates of the unknown parameters by using the 'Params' name-value argument. Otherwise, irf issues an error.

irf issues an error when Mdl is a dimension-varying model, which is a time-varying model containing at least one variable that changes dimension during the sampling period (for example, a state variable drops out of the model).

Tip

If Mdl is fully specified, you cannot estimate confidence bounds. To estimate confidence bounds:

  1. Create a partially specified state-space model template for estimation Mdl.

  2. Estimate the model by using the estimate function and data. Return the estimated parameters estParams and estimated parameter covariance matrix EstParamCov.

  3. Pass the model template for estimation Mdl to irf, and specify the parameter estimates and covariance matrix by using the 'Params' and 'EstParamCov' name-value arguments.

  4. For the irf function, return the appropriate output arguments for lower and upper confidence bounds.

Name-Value Arguments

Specify optional pairs of arguments as Name1=Value1,...,NameN=ValueN, where Name is the argument name and Value is the corresponding value. Name-value arguments must appear after other arguments, but the order of the pairs does not matter.

Before R2021a, use commas to separate each name and value, and enclose Name in quotes.

Example: 'NumPeriods',10,'Cumulative',true specifies a 10-period cumulative IRF starting at time 1, during which irf applies the shock to a state-disturbance variable in the system, and ending at period 10.

IRF Options

collapse all

Number of periods for which irf computes the IRF, specified as a positive integer. Periods in the IRF start at time 1 and end at time NumPeriods.

Example: 'NumPeriods',10 specifies the inclusion of 10 consecutive time points in the IRF starting at time 1, during which irf applies the shock, and ending at time 10.

Data Types: double

Estimates of the unknown parameters in the partially specified state-space model Mdl, specified as a numeric vector.

If Mdl is partially specified (contains unknown parameters specified by NaNs), you must specify Params. The estimate function returns parameter estimates of Mdl in the appropriate form. However, you can supply custom estimates by arranging the elements of Params as follows:

  • If Mdl is an explicitly created model (Mdl.ParamMap is empty []), arrange the elements of Params to correspond to hits of a column-wise search of NaNs in the state-space model coefficient matrices, initial state mean vector, and covariance matrix.

    • If Mdl is time invariant, the order is A, B, C, D, Mean0, and Cov0.

    • If Mdl is time varying, the order is A{1} through A{end}, B{1} through B{end}, C{1} through C{end}, D{1} through D{end}, Mean0, and Cov0.

  • If Mdl is an implicitly created model (Mdl.ParamMap is a function handle), the first input argument of the parameter-to-matrix mapping function determines the order of the elements of Params.

If Mdl is fully specified, irf ignores Params.

Example: Consider the state-space model Mdl with A = B = [NaN 0; 0 NaN] , C = [1; 1], D = 0, and initial state means of 0 with covariance eye(2). Mdl is partially specified and explicitly created. Because the model parameters contain a total of four NaNs, Params must be a 4-by-1 vector, where Params(1) is the estimate of A(1,1), Params(2) is the estimate of A(2,2), Params(3) is the estimate of B(1,1), and Params(4) is the estimate of B(2,2).

Data Types: double

Flag for computing the cumulative IRF, specified as a value in this table.

ValueDescription
trueirf computes the cumulative IRF of all variables over the specified time range.
falseirf computes the standard, period-by-period IRF of all variables over the specified time range.

Example: 'Cumulative',true

Data Types: logical

IRF estimation algorithm, specified as 'repeated-multiplication' or 'eigendecomposition'.

The IRF estimator of time m contains the factor Am. This table describes the supported algorithms to compute the matrix power.

ValueDescription
'repeated-multiplication'irf uses recursive multiplication.
'eigendecomposition'irf attempts to use the spectral decomposition of A to compute the matrix power. Specify this value only when you suspect that the recursive multiplication algorithm might experience numerical issues. For more details, see Algorithms.

Data Types: string | char

Confidence Bound Estimation Options

collapse all

Estimated covariance matrix of unknown parameters in the partially specified state-space model Mdl, specified as a positive semidefinite numeric matrix.

estimate returns the estimated parameter covariance matrix of Mdl in the appropriate form. However, you can supply custom estimates by setting EstParamCov(i,j) to the estimated covariance of the estimated parameters Params(i) and Params(j), regardless of whether Mdl is time invariant or time varying.

If Mdl is fully specified, irf ignores EstParamCov.

By default, irf does not estimate confidence bounds.

Data Types: double

Number of Monte Carlo sample paths (trials) to generate to estimate confidence bounds, specified as a positive integer.

Example: 'NumPaths',5000

Data Types: double

Confidence level for the confidence bounds, specified as a numeric scalar in the interval [0,1].

For each period, randomly drawn confidence intervals cover the true response 100*Confidence% of the time.

The default value is 0.95, which implies that the confidence bounds represent 95% confidence intervals.

Example: Confidence=0.9 specifies 90% confidence intervals.

Data Types: double

Output Arguments

collapse all

IRFs of the measurement variables yt, returned as a NumPeriods-by-k-by-n numeric array.

ResponseY(t,i,j) is the dynamic response of measurement variable j at period t, when a unit shock is applied to state-disturbance variable i during period 1, for t = 1,2,...,NumPeriods, i = 1,2,...,k, and j = 1,2,...,n.

IRFs of the state variables xt, returned as a NumPeriods-by-k-by-m numeric array.

ResponseX(t,i,j) is the dynamic response of state variable j at period t, when a unit shock is applied to state-disturbance variable i during period 1, for t = 1,2,...,NumPeriods, i = 1,2,...,k, and j = 1,2,...,m.

Pointwise lower confidence bounds of the measurement variable IRF, returned as a NumPeriods-by-k-by-n numeric array.

LowerY(t,i,j) is the lower bound of the 100*Confidence% percentile interval on the true dynamic response of measurement variable j at period t, when a unit shock is applied to state-disturbance variable i during period 1.

Pointwise upper confidence bounds of the measurement variable IRF, returned as a NumPeriods-by-k-by-n numeric array.

UpperY(t,i,j) is the upper confidence bound corresponding to the lower confidence bound LowerY(t,i,j).

Pointwise lower confidence bounds of the state variable IRF, returned as a NumPeriods-by-k-by-m numeric array.

LowerX(t,i,j) is the lower bound of the 100*Confidence% percentile interval on the true dynamic response of state variable j at period t, when a unit shock is applied to state-disturbance variable i during period 1.

Pointwise upper confidence bounds of the state variable IRF, returned as a NumPeriods-by-k-by-m numeric array.

UpperX(t,i,j) is the upper confidence bound corresponding to the lower bound LowerX(t,i,j).

More About

collapse all

Impulse Response Function

An impulse response function (IRF) of a state-space model (or dynamic response of the system) measures contemporaneous and future changes in the state and measurement variables when each state-disturbance variable is shocked by a unit impulse at period 1. In other words, the IRF at time t is the derivative of each state and measurement variable at time t with respect to a state-disturbance variable at time 1, for each t ≥ 1.

Consider the time-invariant state-space model

xt=Axt1+Butyt=Cxt+Dεt,

and consider an unanticipated unit shock at period 1, applied to state-disturbance variable j uj,t.

The r-step-ahead response of the state variables xt to the shock is

ψxj(r)=Arbj,

where r > 0 and bj is column j of the state-disturbance-loading matrix B.

The r-step-ahead response of the measurement variables yt to the shock is

ψyj(r)=CArbj.

IRFs depend on the time interval over which they are computed. However, the IRF of a time-invariant state-space model is time homogeneous, which means that the IRF does not depend on the time at which the shock is applied. Time-varying IRFs, which are the IRFs of a time-varying but dimension-invariant system, have the form

ψxj(r)=ArA2A1b1,jψyj(r)=CrArA2A1b1,j,

where b1,j is column j of B1, the period 1 state-disturbance-loading matrix. Time-varying IRFs depend on the time at which the shock is applied. irf always applies the shock at period 1.

IRFs are independent of the initial state distribution.

Algorithms

  • If you specify 'eigendecomposition' for the 'Method' name-value pair argument, irf attempts to diagonalize the state-transition matrix A by using the spectral decomposition. irf resorts to recursive multiplication instead under at least one of these circumstances:

    • An eigenvalue is complex.

    • The rank of the matrix of eigenvectors is less than the number of states

    • Mdl is time varying.

  • irf uses Monte Carlo simulation to compute confidence intervals.

    1. irf randomly draws NumPaths variates from the asymptotic sampling distribution of the unknown parameters in Mdl, which is Np(Params,EstParamCov), where p is the number of unknown parameters.

    2. For each randomly drawn parameter set j, irf:

      1. Creates a state-space model that is equal to Mdl, but substitutes in parameter set j

      2. Computes the random IRF of the resulting model ψj(t), where t = 1 through NumPaths

    3. For each time t, the lower bound of the confidence interval is the (1 – c)/2 quantile of the simulated IRF at period t ψ(t), where c = Confidence. Similarly, the upper bound of the confidence interval at time t is the (1 – c)/2 upper quantile of ψ(t).

Version History

Introduced in R2020b