*
* BOOTVAR.PRG
* Bootstrapping example for VAR
*
compute lags=4			   ;*Number of lags
compute nvar=2			   ;*Number of variables
compute nstep=10		   ;*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) responsessq(nkeep)
declare rect[series] impulses(nvar,nvar) impulsesc(nvar,nvar) impulsessq(nvar,nvar)

*
infobox(action=define,progress,lower=1,upper=nkeep) "Bootstrap Simulations"
do draws=1,nkeep
   *
   *
   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(noprint,noftests,cvout=v)
   impulse(noprint,model=bootvar,cv=v,results=impulses,steps=nstep)




    *TAKING SQUARES OF THE IMPULSES
    do i=1, nvar
   do j=1, nvar
      set impulsessq(i,j) =impulses(i,j)^2
   end do j
   end do i

   *
   dim responses(draws)(nvar*nvar,nstep)
   dim responsessq(draws)(nvar*nvar,nstep)

  ewise responses(draws)(i,j)=impulses((i-1)/nvar+1,%clock(i,nvar))(j)
  ewise responsessq(draws)(i,j)=impulsessq((i-1)/nvar+1,%clock(i,nvar))(j)

   infobox(current=draws)
end do draws
infobox(action=remove)

*


*1. reading the squarred response functions, which are squared inside the bootstrap loop
dec vect[series] upper(nvar) lower(nvar) resp(nvar)
dec rect[series] respsq(nvar,nvar)
do i=1,nvar
   smpl 1 nkeep
   do j=1,nvar
      clear  resp(j)
      do k=1,nstep
         set work   1 nkeep = responsessq(t)((i-1)*nvar+j,k)
         compute resp(j)(k)=%avg(work)
      end do k
   end do j
   smpl 1 nstep
   do j=1,nvar
      set respsq(i,j) = resp(j)
   end do j
end do i

* 2. reading the response function
dec rect[series] resps(nvar,nvar)
do i=1,nvar
   smpl 1 nkeep
   do j=1,nvar
      clear  resp(j)
      do k=1,nstep
         set work   1 nkeep = responses(t)((i-1)*nvar+j,k)
         compute resp(j)(k)=%avg(work)
      end do k
   end do j
   smpl 1 nstep
   do j=1,nvar
      set resps(i,j) = resp(j)
   end do j
end do i

*3. the difference between the response functions, which are squared outside the bootstrap loop and those, which are squared inside the loop should be zero
set con = resps(2,2)^2-respsq(2,2)
print / con