*
* BOOTVAR.PRG
* Bootstrapping example for VAR
*
compute lags=4			   ;*Number of lags
compute nvar=2			   ;*Number of variables
compute nstep=16		   ;*Number of response steps
compute nkeep=2500		;*Number of draws
*
open data haversample.rat
calendar(q) 1959
data(format=rats) 1959:1 2006:4 fm1 gdph
*
log gdph
log fm1
******************************************************************
*
* Set up the system
*
dec vect[series] udraws(nvar) resids(nvar) resample(nvar)
dec vect[equation] eqsample(nvar) eqbase(nvar)
*
system(model=varmodel)
variables gdph fm1
lags 1 to lags
det constant
end(system)
*
estimate(resids=resids)
compute rstart=%regstart(),rend=%regend()
*
* Set up the parallel system for the resampled data
*
system(model=bootvar)
variables resample
lags 1 to lags
det constant
end(system)
*
do i=1,nvar
   set resample(i) = %modeldepvars(varmodel)(i){0}
end do i
*
declare vect[rect] responses(nkeep)
declare rect[series] impulses(nvar,nvar)
*
infobox(action=define,progress,lower=1,upper=nkeep) "Bootstrap Simulations"
do draws=1,nkeep
   *
   * Draw a bootstrapped sample
   *
   boot entries rstart rend
   do i=1,nvar
      set udraws(i) = resids(i)(entries(t))
   end do i
   forecast(paths,model=varmodel,from=rstart,to=rend,results=resample)
   # udraws
   *
   * Estimate the model with resampled data
   *
   estimate(noprint,noftests,cvout=v)
   impulse(noprint,model=bootvar,cv=v,results=impulses,steps=nstep)
   *
   * Store the impulse responses
   *
   dim responses(draws)(nvar*nvar,nstep)
   ewise responses(draws)(i,j)=impulses((i-1)/nvar+1,%clock(i,nvar))(j)
   infobox(current=draws)
end do draws
infobox(action=remove)

dec vect[strings] xlabel(nvar) ylabel(nvar)
dec vect[integer] depvars
compute depvars=%modeldepvars(varmodel)
do i=1,nvar
   compute ll=%l(depvars(i))
   compute xlabel(i)=ll
   compute ylabel(i)=ll
end do i

grparm(bold) hlabel 18 matrixlabels 14
grparm  axislabel 24
spgraph(header="Impulse responses",xpos=both,xlab=xlabel, $
        ylab=ylabel,vlab="Responses of",vfields=nvar,hfields=nvar)

*
* Because we want a common scale for all responses of a single variable, we need
* to do all the calculations for a full row of graphs first. Note, by the way,
* that this graph transposes rows and columns from the arrangement used in the
* original montevar.
*

dec vect[series] upper(nvar) lower(nvar) resp(nvar)
do i=1,nvar
   compute minlower=maxupper=0.0
   smpl 1 nkeep
   do j=1,nvar
      clear lower(j) upper(j) resp(j)
      do k=1,nstep
         set work   1 nkeep = responses(t)((i-1)*nvar+j,k)
         compute frac=%fractiles(work,||.16,.84||)
         compute lower(j)(k)=frac(1)
         compute upper(j)(k)=frac(2)
         compute resp(j)(k)=%avg(work)
      end do k
      compute maxupper=%max(maxupper,%maxvalue(upper(j)))
      compute minlower=%min(minlower,%minvalue(lower(j)))
   end do j
*
   smpl 1 nstep
   do j=1,nvar
      graph(ticks,min=minlower,max=maxupper,number=0) 3 j i
      # resp(j)
      # upper(j) / 2
      # lower(j) / 2
   end do j
end do i
*
*
spgraph(done)