Stata The Stata listserver
[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

st: Problem with CIRF & COIRF in var


From   "Cavallo, Alexander" <acavallo@lexecon.com>
To   "'Statalist (statalist@hsphsun2.harvard.edu)'" <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/



© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index