/****************** Time Series Standard Errors **********************
Author: Thierry.Roncalli Thierry.Roncalli@montesquieu.u-bordeaux.fr
Date: 25 Oct 1996
This code is provided without guarantee for public non-commercial use.

Purpose: compute the standard errors for forecast error
impulse, orthogonalized impulse and forecast error variance
decomposition.

A possible answer is to use the Monte Carlo method:
For example, we can simulate different paths of the estimates,
compute the impulse function for each path and then compute
the standard error.
Another possibility is to use the bootstrap method of Stoffer and Wall.
We simulate different paths of the data set with the bootstrap
method. For each path, we compute the estimates. Then, we compute
the standard error.
The two following examples reproduces the figures 3.4, 3.5, 3.6
and 3.7 of LUTKEPOHL [1991]. The example Impulse1.prg implement
the first method (5000 simulations take 75 seconds with a pentium 90).
The example Impulse2.prg consider the bootstrap
method of Stoffer and Wall (500 simulations take 210 seconds
with a pentium 90). We use the arma_to_ssm procedure
to convert the VAR model into a SSM model. Then, we simulate a
path with the bootstrap_SSM procedure. Note that
we can use these two methods to estimate the standard error
of the impulse function for a vector arma model or a state
space model.

Requires the TSM package.
*********************************************************************/

@Impulse1.prg --->@

/*
**  LUTKEPOHL [1991], Introduction to Multiple Time Series Analysis,
**  Springer-Verlag, Berlin-Heidelberg
**
**  Impulse Response Analysis
**  See LUTKEPOHL, page 107
**
*/

new;
library tsm,optmum,pgraph;
TSMset;

load y[92,3] = lutkepoh.asc;

y = ln(y[1:76,.]);

INVEST = y[.,1];
dinv = INVEST-lag1(INVEST);
INCOME = y[.,2];
dinc = INCOME-lag1(INCOME);
CONSUM = y[.,3];
dcons = CONSUM-lag1(CONSUM);
data = dinv~dinc~dcons;

_print = 0;
{theta,stderr,Mcov,LOGL} = varx_LS(data,1,2);  /* Constant and 2 lags */

/* AR part of the VARX model */

theta = theta[1:18];
Mcov = Mcov[1:18,1:18];

Nr = 8; /* Number of responses */

beta = theta;
call arma_impulse(beta,2,0,Nr);

rep1 = zeros(1+Nr,1);  /* Responses of Consumption to an
                          impulse in Income */
rep2 = zeros(1+Nr,1);  /* Responses of Investment to an
                          impulse in Consumption */
rep3 = zeros(1+Nr,1);  /* Accumulated responses of Consumption to an
                          impulse in Income */
rep4 = zeros(1+Nr,1);  /* Accumulated responses of Investment to an
                          impulse in Consumption */


j = 0;
do until j > Nr;
  z = varget("IMPULSE"$+ftos(j,"%lf",1,0));
  rep1[1+j,1] = z[3,2];
  rep2[1+j,1] = z[1,3];
  z = varget("_IMPULSE"$+ftos(j,"%lf",1,0));
  rep3[1+j,1] = z[3,2];
  rep4[1+j,1] = z[1,3];
  j = j + 1;
endo;

Ns = 5000;

RND_theta = rndn2(theta,Mcov,Ns);

rep1s = zeros(1+Nr,Ns);
rep2s = zeros(1+Nr,Ns);
rep3s = zeros(1+Nr,Ns);
rep4s = zeros(1+Nr,Ns);


i = 1;
do until i > Ns;

  beta = RND_theta[i,.]';
  call arma_impulse(beta,2,0,Nr);

  j = 0;
  do until j > Nr;
    z =  varget("IMPULSE"$+ftos(j,"%lf",1,0));
    rep1s[1+j,i] = z[3,2];
    rep2s[1+j,i] = z[1,3];
    z =  varget("_IMPULSE"$+ftos(j,"%lf",1,0));
    rep3s[1+j,i] = z[3,2];
    rep4s[1+j,i] = z[1,3];
    j = j + 1;
  endo;

  i = i + 1;
endo;

rep1s = rep1s';
stderr1 = stdc(rep1s);
rep2s = rep2s';
stderr2 = stdc(rep2s);
rep3s = rep3s';
stderr3 = stdc(rep3s);
rep4s = rep4s';
stderr4 = stdc(rep4s);

t = seqa(0,1,Nr+1);

graphset;
  _pnum = 2; _pdate = ""; _pltype = 1|6|1; _pframe = {0,0}; _pcross = 1;
  bound = rep1 + (-1.96~0~1.96).*stderr1;
  xy(seqa(0,1,Nr+1),bound);

graphset;
  _pnum = 2; _pdate = ""; _pltype = 1|6|1; _pframe = {0,0}; _pcross = 1;
  bound = rep2 + (-1.96~0~1.96).*stderr2;
  xy(seqa(0,1,Nr+1),bound);

graphset;
  _pnum = 2; _pdate = ""; _pltype = 1|6|1; _pframe = {0,0}; _pcross = 1;
  bound = rep3 + (-1.96~0~1.96).*stderr3;
  xy(seqa(0,1,Nr+1),bound);

graphset;
  _pnum = 2; _pdate = ""; _pltype = 1|6|1; _pframe = {0,0}; _pcross = 1;
  bound = rep4 + (-1.96~0~1.96).*stderr4;
  xy(seqa(0,1,Nr+1),bound);


proc (1) = rndn2(mu,SIGMA,Ns);
  local d,u,Pchol,y,i,ui;
  d = rows(mu);
  u = rndn(d,Ns);
  Pchol = chol(SIGMA)';
  y = zeros(d,Ns);
  i = 1;
  do until i > Ns;
    ui = u[.,i];
    y[.,i] = Pchol*ui;
    i=i+1;
  endo;
  y = mu + y;
  y = y';
  retp(y);
endp;



@Impulse2.prg --->@

/*
**  LUTKEPOHL [1991], Introduction to Multiple Time Series Analysis,
**  Springer-Verlag, Berlin-Heidelberg
**
**  Impulse Response Analysis
**  See LUTKEPOHL, page 107
**
*/

new;
library tsm,optmum,pgraph;
TSMset;

load y[92,3] = lutkepoh.asc;

y = ln(y[1:76,.]);

INVEST = y[.,1];
dinv = INVEST-lag1(INVEST);
INCOME = y[.,2];
dinc = INCOME-lag1(INCOME);
CONSUM = y[.,3];
dcons = CONSUM-lag1(CONSUM);
{data,retcode} = Missing(dinv~dinc~dcons,0);
Nobs = rows(data);

_print = 0;
{theta,stderr,Mcov,LOGL} = varx_LS(data,1,2);  /* Constant and 2 lags */

/* AR part of the VARX model */

Nr = 8; /* Number of responses */

beta = theta[1:18];
const = theta[19:21];
SIGMA = _varx_SIGMA;

call arma_impulse(beta,2,0,Nr);

rep1 = zeros(1+Nr,1);  /* Responses of Consumption to an
                          impulse in Income */
rep2 = zeros(1+Nr,1);  /* Responses of Investment to an
                          impulse in Consumption */
rep3 = zeros(1+Nr,1);  /* Accumulated responses of Consumption to an
                          impulse in Income */
rep4 = zeros(1+Nr,1);  /* Accumulated responses of Investment to an
                          impulse in Consumption */


j = 0;
do until j > Nr;
  z = varget("IMPULSE"$+ftos(j,"%lf",1,0));
  rep1[1+j,1] = z[3,2];
  rep2[1+j,1] = z[1,3];
  z = varget("_IMPULSE"$+ftos(j,"%lf",1,0));
  rep3[1+j,1] = z[3,2];
  rep4[1+j,1] = z[1,3];
  j = j + 1;
endo;

Ns = 500;

rep1s = zeros(1+Nr,Ns);
rep2s = zeros(1+Nr,Ns);
rep3s = zeros(1+Nr,Ns);
rep4s = zeros(1+Nr,Ns);


i = 1;
do until i > Ns;

  {Z,d,H,T,c,R,Q} = arma_to_SSM(beta,2,0,SIGMA);
  c[1:3] = const;
  call SSM_build(Z,d,H,T,c,R,Q,0);
  a0 = data[1,.]'|data[2,.]'|zeros(3,1);
  P0 = zeros(9,9);
  call KFiltering(data[3:Nobs,.],a0,P0);
  {ys,as} = bootstrap_SSM(a0);
  datas = data[1,.]|data[2,.]|ys;
  {RND_theta,stderr,Mcov,LOGL} = varx_LS(datas,1,2);
  betas = RND_theta[1:18];
  call arma_impulse(betas,2,0,Nr);

  j = 0;
  do until j > Nr;
    z =  varget("IMPULSE"$+ftos(j,"%lf",1,0));
    rep1s[1+j,i] = z[3,2];
    rep2s[1+j,i] = z[1,3];
    z =  varget("_IMPULSE"$+ftos(j,"%lf",1,0));
    rep3s[1+j,i] = z[3,2];
    rep4s[1+j,i] = z[1,3];
    j = j + 1;
  endo;

  i = i + 1;
endo;

rep1s = rep1s';
stderr1 = stdc(rep1s);
rep2s = rep2s';
stderr2 = stdc(rep2s);
rep3s = rep3s';
stderr3 = stdc(rep3s);
rep4s = rep4s';
stderr4 = stdc(rep4s);

t = seqa(0,1,Nr+1);

graphset;
  _pnum = 2; _pdate = ""; _pltype = 1|6|1; _pframe = {0,0}; _pcross = 1;
  bound = rep1 + (-1.96~0~1.96).*stderr1;
  xy(seqa(0,1,Nr+1),bound);

graphset;
  _pnum = 2; _pdate = ""; _pltype = 1|6|1; _pframe = {0,0}; _pcross = 1;
  bound = rep2 + (-1.96~0~1.96).*stderr2;
  xy(seqa(0,1,Nr+1),bound);

graphset;
  _pnum = 2; _pdate = ""; _pltype = 1|6|1; _pframe = {0,0}; _pcross = 1;
  bound = rep3 + (-1.96~0~1.96).*stderr3;
  xy(seqa(0,1,Nr+1),bound);

graphset;
  _pnum = 2; _pdate = ""; _pltype = 1|6|1; _pframe = {0,0}; _pcross = 1;
  bound = rep4 + (-1.96~0~1.96).*stderr4;
  xy(seqa(0,1,Nr+1),bound);