# st: Problem with CIRF & COIRF in var

 From "Cavallo, Alexander" To "'Statalist (statalist@hsphsun2.harvard.edu)'" Subject st: Problem with CIRF & COIRF in var Date Tue, 9 Nov 2004 15:26:49 -0600

```I am trying to replicate the CIRF and COIRF calculations for infinite
horizon in the Lutkepohl 1993 book.  I am matching on the calculation of the
infinite CIRF and COIRF but not on the asymptotic standard errors.  See pp
97-99 in "Introduction to Multiple Time Series".  Estimates of CIRF and
COIRF should match pp 106 and 111.

I have tried to hack the VARIRF commands but could not figure out how or
where to modify it.  I would appreciate any corrections to this code.

Thanks!!

--Alex Cavallo
Lexecon, Inc.

program define varirf_inf
mat a=e(b)
mat V_a=e(V)
mat sigma_u=e(Sigma)
local k=e(neqs)
local p=e(mlag)
local T=e(N)

* drop constants
mat R0=(I(`k'*`p'), J(`k'*`p',`k'*`k'*`p'+`k'-`k'*`p',0))
local t=1
while `t'<`k'-1 {
local pre=`t'*(`k'*`p'+1)
local post=`k'*`k'*`p'+`k'-`k'*`p'-`pre'
local t=`t'+1
mat R=(J(`k'*`p',`pre',0), I(`k'*`p'), J(`k'*`p',`post',0))
mat R0=R0\ R
}
mat R=(J(`k'*`p',`k'*`k'*`p'+`k'-`k'*`p'-1,0), I(`k'*`p'), J(`k'*`p',1,0))
mat R0=R0\ R
mat a=R0*a'
mat V_a=R0*V_a*R0'

* reshape var coefs into Lutkepohl format
mat R=J(`k'*`k'*`p',`k'*`k'*`p',0)
local lag=0
local ctr=0
while `lag'<`p' {
local lag=`lag'+1
local s=0
while `s'<`k' {
local s=`s'+1
local t=0
while `t'<`k' {
local t=`t'+1
local ctr=`ctr'+1
local z=`lag'+(`s'-1)*`p'+(`t'-1)*`k'*`p'
mat R[`ctr',`z']=1
}
}
}
mat alpha=R*a
mat sigma_a=R*V_a*R'

* make J matrix
mat j=J(`k',`k'*`p',0)
local t=0
while `t'<`k' {
local t=`t'+1
mat j[`t',`t']=1
}

* make elimination matrix
local lr=`k'*(`k'+1)/2
local lc=`k'*`k'
local t=`k'-1
local s=0
local fact=`k'+1
mat L=I(`k'),J(`k',`lc'-`k',0)
while `t'>1 {
local s=`s'+1
mat l=J(`t',`s'*`fact',0),I(`t'),J(`t',`lc'-`t'-`s'*`fact',0)
mat L=L\l
local t=`t'-1
}
local s=`s'+1
mat l=J(`t',`s'*`fact',0),I(`t')
mat L=L\l

* make commutation matrix
mat K=J(`k'*`k',`k'*`k',0)
local i=0
local t=0
while `i'<`k' {
local i=`i'+1
local j=0
while `j'<`k' {
local j=`j'+1
local t=`t'+1
local s=(`j'-1)*`k'+`i'
di `t' `s'
mat K[`t',`s']=1
}
}

* make duplication matrix
matrix D = J(`k'*`k',round(.5*(`k'+1)*`k',1),0)

local cin  1
forvalues j=1(1)`k' {
forvalues i=1(1)`k' {
local r = `i' + (`j'-1)*`k'
if `i'>=`j' {
local c = `cin'
local cin = `cin' + 1
}
else {
local q = round( .5*( `k'*(`k'+1) /*
*/ - (`k'-`i'+1)*(`k'-`i'+2) ),1)
local c = `q' + `j'-(`i'-1)
}
matrix D[`r',`c']=1
}

}

* var matrix of error estimates
mat sigma_u=e(Sigma)
mat isu=inv(sigma_u)
mat sigma_s=2*inv(D'*(isu#isu)*D)
mat li sigma_s

* orthogonalize errors
mat P=cholesky(sigma_u)

* make Ai matrices
local l=0
local z=0
while `l'<`p' {
local l=`l'+1
mat A`l'=J(`k',`k',0)
local s=0
while `s'<`k' {
local s=`s'+1
local t=0
while `t'<`k' {
local t=`t'+1
local z=`z'+1
mat A`l'[`t',`s']=alpha[`z',1]
}
}
}

* compute matrices for infinte cumulative IRF
mat CIRFinf=I(`k')
local l=0
while `l'<`p' {
local l=`l'+1
mat CIRFinf=CIRFinf-A`l'
}
mat CIRFinf=inv(CIRFinf)
mat COIRFinf=CIRFinf*P

mat Finf=CIRFinf'
local l=1
while `l'<`p' {
local l=`l'+1
mat Finf=Finf,CIRFinf'
}
mat Finf=Finf#CIRFinf
mat Binf=(P'#I(`k'))*Finf
mat H=L'*inv( L*(I(`k'*`k')+K)*(P#I(`k'))*L')
mat Bbarinf=(I(`k')#CIRFinf)*H

mat V_COIRFinf=(Binf*sigma_a*Binf' + Bbarinf*sigma_s*Bbarinf')
mat se=vecdiag(V_COIRFinf)
local t=0
while `t'<`k'*`k' {
local t=`t'+1
mat se[1,`t']=sqrt(se[1,`t'])
}
mat se=se'
mat SE_COIRFinf=J(`k',`k',0)
local s=0
local z=0
while `s'<`k' {
local s=`s'+1
local t=0
while `t'<`k' {
local t=`t'+1
local z=`z'+1
mat SE_COIRFinf[`t',`s']=se[`z',1]
}
}
mat V_CIRFinf=(Finf*sigma_a*Finf')
mat se=vecdiag(V_CIRFinf)
local t=0
while `t'<`k'*`k' {
local t=`t'+1
mat se[1,`t']=sqrt(se[1,`t'])
}
mat se=se'
mat SE_CIRFinf=J(`k',`k',0)
local s=0
local z=0
while `s'<`k' {
local s=`s'+1
local t=0
while `t'<`k' {
local t=`t'+1
local z=`z'+1
mat SE_CIRFinf[`t',`s']=se[`z',1]
}
}

mat li CIRFinf
mat li SE_CIRFinf

mat li COIRFinf
mat li SE_COIRFinf
end

clear
set more off
capture log close
log using lutk2.log, replace text

webuse lutkepohl
tsset
keep if qtr<=q(1978q4)
var dlinvestment dlincome dlconsumption, lags(1/2) lutstats dfk
varirf create lutk, step(100) set(lutk) replace
varirf table cirf

varirf_inf
*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/
```