Kilian 1998 (RES)
Posted: Wed Jul 27, 2011 10:52 amHallo,
I tried to replicate Kilian (1998), Small-Sample Conficence Intervals for Impulse Response Functions, RES (http://www.ssc.wisc.edu/~bhansen/718/Kilian1998.pdf). This is a procedure that takes a VAR model and computes confidence intervals for the IRFs based on the bootstrap after bootstrap approach explained in the paper (plus an auxiliary function):
The second file uses this procedure and replicates the MC-simulation shown in the paper for one specific parameterization of the DGP:
I suspect that there is some bug in the code because the intervals obtained by Kilian's method do almost never include the IRFs obtained from the true model. Thus, if anybody is also interested in this approach he/she might have a look at the code to try to find out whether there is something coded totally wrong or the strange results are just due to a syntax error that I overlooked.
I tried to replicate Kilian (1998), Small-Sample Conficence Intervals for Impulse Response Functions, RES (http://www.ssc.wisc.edu/~bhansen/718/Kilian1998.pdf). This is a procedure that takes a VAR model and computes confidence intervals for the IRFs based on the bootstrap after bootstrap approach explained in the paper (plus an auxiliary function):
- Code: Select all
function %%vecvecavg array
type vec[rea] %%vecvecavg
type vec[vec[rea]] array
*
local integer n k i j
local vec temp
*
comp n = %rows(array)
comp k = %rows(array(1))
dim temp(n) %%vecvecavg(k)
do i=1,k
ewise temp(j) = array(j)(i)
comp %%vecvecavg(i) = %avg(temp)
end do i
end function
************************************************
procedure KILIAN98 start end
type integer start end
*
option model model
option integer ndraws 1000
option integer niters 2000
option real signif .95
option integer steps 8
option switch graph 0
option switch display 0
*
local integer nvar i j k draw time
local vec[ser] origser resids udraws artifser upper lower median
local rec[ser] girf impulses
local vec betaols bias betatilde betatildestar
local vec[com] cv
local ser[int] entries
local vec[vec[rea]] betastar
local vec[rec[ser]] bgirf
local vec work frac
local real maxupper minlower adjfactor
local model origmodel
*
comp nvar = %modelsize(model)
dim udraws(nvar) betastar(ndraws) bgirf(niters) upper(nvar) lower(nvar) median(nvar) girf(nvar,nvar)
do draw=1,niters
dim bgirf(draw)(nvar,nvar)
end do i
*
* Save original data
dim origser(nvar)
do i=1,nvar
set origser(i) start end = ([series] %modeldepvars(model)(i))
end do i
*
* Estimate model and save IRFs based on OLS estimates
estimate(noprint,model=model,resid=resids) start end
comp origmodel = model
comp betaols = %BETASYS
impulse(noprint,model=model,results=girf,steps=steps)
*
infobox(action=define,progress,lower=1,upper=ndraws) "Estimating Bias"
do draw = 1,ndraws
infobox(current=draw)
*
* Restore original data and model
do i=1,nvar
set ([series] %modeldepvars(model)(i)) start end = origser(i)(t)
end do i
comp %modelsetcoeffs(model,betaols)
*
* Draw a bootstrapped sample
boot entries start end start end
do i=1,nvar
set udraws(i) start end = resids(i)(entries(t))
end do i
*
forecast(paths,model=model,from=start,to=end,nostatic,results=artifser,noprint)
# udraws
*
do i=1,nvar
set ([series] %modeldepvars(model)(i)) start end = artifser(i)(t)
end do i
*
* Estimate VAR based on resampled data
estimate(noprint,model=model) start end
comp betastar(draw) = %BETASYS
end do draw
infobox(action=remove)
*
* Compute estimate of bias
comp bias = %%vecvecavg(betastar)-betaols
*
* Check stationarity and correct for bias
comp %modelsetcoeffs(model,betaols)
eigen(cvalues=cv) %modelcompanion(model)
if %cabs(cv(1))>=1.
comp betatilde = betaols
else {
comp adjfactor=1., betatilde = betaols-adjfactor*bias
comp %modelsetcoeffs(model,betatilde)
eigen(cvalues=cv) %modelcompanion(model)
while %cabs(cv(1))>=1. {
comp betatilde = betaols-adjfactor*bias
comp %modelsetcoeffs(model,betatilde)
eigen(cvalues=cv) %modelcompanion(model)
comp adjfactor=adjfactor-0.01
}
}
*
* Using betatilde to perform the second bootstrap
comp %modelsetcoeffs(model,betatilde)
forecast(nopaths,model=model,from=start,to=end,nostatic,results=artifser,noprint)
do i=1,nvar
set resids(i) start end = origser(i)(t)-artifser(i)(t)
end do i
*
dim betastar(niters)
infobox(action=define,progress,lower=1,upper=niters) "Running Second Bootstrap"
do draw = 1,niters
infobox(current=draw)
*
* Restore original data and model (with betatilde)
do i=1,nvar
set ([series] %modeldepvars(model)(i)) start end = origser(i)(t)
end do i
comp %modelsetcoeffs(model,betatilde)
*
* Draw a bootstrapped sample
boot entries start end start end
do i=1,nvar
set udraws(i) start end = resids(i)(entries(t))
end do i
*
forecast(paths,model=model,from=start,to=end,nostatic,results=artifser,noprint)
# udraws
*
do i=1,nvar
set ([series] %modeldepvars(model)(i)) start end = artifser(i)(t)
end do i
*
* Estimate VAR based on resampled data
estimate(noprint,model=model) start end
comp betastar(draw) = %BETASYS
*
* Check stationarity and correct for bias
comp %modelsetcoeffs(model,betatilde)
eigen(cvalues=cv) %modelcompanion(model)
if %cabs(cv(1))>=1.
comp betatildestar = betatilde
else {
comp adjfactor=1., betatildestar = betatilde-adjfactor*bias
comp %modelsetcoeffs(model,betatildestar)
eigen(cvalues=cv) %modelcompanion(model)
while %cabs(cv(1))>=1. {
comp betatildestar = betatilde-adjfactor*bias
comp %modelsetcoeffs(model,betatildestar)
eigen(cvalues=cv) %modelcompanion(model)
comp adjfactor=adjfactor-0.01
}
}
impulse(noprint,model=model,results=impulses,steps=steps)
*
do i=1,nvar
do j=1,nvar
set bgirf(draw)(i,j) 1 steps = impulses(i,j)(t)
end do j
end do i
end do draw
infobox(action=remove)
*
* Restore original data
do i=1,nvar
set ([series] %modeldepvars(model)(i)) start end = origser(i)(t)
end do i
*
* Compute fractiles of IRFs
decl rect[ser] %UPPER(nvar,nvar) %LOWER(nvar,nvar)
dim work(niters)
do i=1,nvar
do j=1,nvar
do time=1,steps
ewise work(k) = bgirf(k)(i,j)(time)
comp frac=%fractiles(work,||(1-signif)/2.,.5,(1+signif)/2.||)
set upper(j) time time = frac(3)
set median(j) time time = frac(2)
set lower(j) time time = frac(1)
end do time
set %UPPER(i,j) 1 steps = upper(j)(t)
set %LOWER(i,j) 1 steps = lower(j)(t)
end do j
end do i
*
* Plot results
if graph {
spgraph(hfi=nvar,vfi=nvar)
do i=1,nvar
comp minlower=maxupper=.0
do j=1,nvar
do time=1,steps
ewise work(k) = bgirf(k)(i,j)(time)
comp frac=%fractiles(work,||(1-signif)/2.,.5,(1+signif)/2.||)
set upper(j) time time = frac(3)
set median(j) time time = frac(2)
set lower(j) time time = frac(1)
end do time
comp maxupper=%max(maxupper,%maxvalue(upper(j)))
comp minlower=%min(minlower,%minvalue(lower(j)))
end do j
*
do j=1,nvar
graph(ticks,min=minlower,max=maxupper,number=0) 4 i j
# girf(i,j) 1 steps;
# median(j) 1 steps 3
# upper(j) 1 steps 2
# lower(j) 1 steps 2
end do j
end do i
spgraph(done)
}
*
if display {
disp; disp; disp " Variable OLS/IV Bootstrap Bias "
do i=1,nvar
disp; disp %l(%eqndepvar(%modeleqn(model,i)))
do j=1,%rows(%modelgetcoeffs(origmodel))
disp j @+0 "." @5 %l(%eqntable(%modeleqn(model,i))(1,j)) @+0 "{" @+0 %eqntable(%modeleqn(model,i))(2,j) @+0 "}" @30 "" $
@@.2 *.##### betaols((i-1)*%rows(%modelgetcoeffs(origmodel))+j) $
@@.10 *.##### betatilde((i-1)*%rows(%modelgetcoeffs(origmodel))+j) $
@@.10 *.##### bias((i-1)*%rows(%modelgetcoeffs(origmodel))+j)
end do j
end do i
disp; disp "Number of replications for bias estimation:" *. ndraws
disp "Number of replications for 2nd bootstrap:" *. niters
}
*
end procedure KILIAN98
The second file uses this procedure and replicates the MC-simulation shown in the paper for one specific parameterization of the DGP:
- Code: Select all
all 2000
comp B = .5
comp sigma = ||1, 0.3|0.3, 1||
comp F = %decomp(sigma)
decl ser[vec] u
gset u = %ranmvnormal(F)
set y1 1 1 = 0
set y2 1 1 = 0
do time=2,2000
set y1 time time = B*y1{1}+u(time)(1)
set y2 time time = 0.5*y1{1}+0.5*y2{1}+u(time)(2)
end do time
graph 2; # y1; # y2
* Set up true model and compute true IRFs
system(model=true)
var y1 y2
det constant
lags 1 to 1
end(system)
estimate
comp %modelsetcoeffs(true,||B,.5|.0,.5|.0,.0||)
impulse(print,model=true,results=trueirf,steps=12)
source(echo) kilian98.src
@kilian98(model=true,ndraws=1000,steps=12,niters=1000,graph,display) 151 2000
* Do MC-Simulation
comp nruns=200
decl vec[rec[ser]] montiupper(nruns) montilower(nruns)
infobox(action=define,progress,lower=1,upper=200)
do run=1,nruns
infobox(current=run)
* Generate data
gset u = %ranmvnormal(F)
set y1 1 1 = 0
set y2 1 1 = 0
do time=2,200
set y1 time time = B*y1{1}+u(time)(1)
set y2 time time = 0.5*y1{1}+0.5*y2{1}+u(time)(2)
end do time
*
* Estimate model
system(model=var)
var y1 y2
det constant
lags 1 to 1
end(system)
@kilian98(model=var,ndraws=500,steps=12,niters=500,nograph,nodisplay) 151 200
dim montiupper(run)(2,2) montilower(run)(2,2)
do i=1,2
do j=1,2
set montiupper(run)(i,j) 1 12 = %UPPER(i,j)(t)
set montilower(run)(i,j) 1 12 = %LOWER(i,j)(t)
end do j
end do i
end do run
infobox(action=remove)
* Get coverage rates
decl rec[ser] coverage(2,2)
do i=1,2
do j=1,2
do step=1,8
comp count=0
do run=1,200
if trueirf(i,j)(step)>=montilower(run)(i,j)(step).and.trueirf(i,j)(step)<=montiupper(run)(i,j)(step)
comp count=count+1
end do run
set coverage(i,j) step step = count/200.*100
end do step
end do j
end do i
graph(header="Coverages for coefficients of A") 4;
# coverage(1,1)
# coverage(1,2)
# coverage(2,1)
# coverage(2,2)
I suspect that there is some bug in the code because the intervals obtained by Kilian's method do almost never include the IRFs obtained from the true model. Thus, if anybody is also interested in this approach he/she might have a look at the code to try to find out whether there is something coded totally wrong or the strange results are just due to a syntax error that I overlooked.