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
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)