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/
