OPEN DATA "C:\Users\n00605626\Desktop\contractionary effect of exchange rate\Gil's data\Austral_2010.xls"
CALENDAR(Q) 1970:01
ALL 2009:04
DATA(FORMAT=XLS,ORG=COLUMNS) 1970:01 2010:04 CAR KAR IP RCP NX RX FGDP FI MS CPI

set CAR = CAR*100

set CPI = log(CPI)
*
system(model=varmodel)
variables  CAR CPI IP RX MS
lags 1 to 4
det constant FGDP{1 to 2} FI{1 to 2}
end(system)
estimate(outsigma=vmat)

*
source(echo) uhligfuncs.src
dec vect[strings] vl(5)
compute vl=||'Current account','consumer price','output', 'exchange rate','money supply'||

*
* n1 is the number of draws from the posterior of the VAR
* n2 is the number of draws from the unit sphere for each draw for the VAR
* nvar is the number of variables
* nstep is the number of IRF steps to compute
* KMAX is the "K" value for the number of steps constrained
*
compute n1=200
compute n2=200
compute nkeep=1000
compute nvar=5
compute nstep=60
compute KMIN=1
compute KMAX=5
*
* This is the standard setup for MC integration of an OLS VAR
*
compute sxx    =%decomp(%xx)
compute svt    =%decomp(inv(%nobs*%sigma))
compute betaols=%modelgetcoeffs(varmodel)
compute ncoef  =%rows(sxx)
compute wishdof=%nobs-ncoef
*
* Most draws are going to get rejected. We allow for up to nkeep good
* ones. The variable accept will count the number of accepted draws.
* GOODRESP will be a RECT(nsteps,nvar) at each accepted draw.
*
declare vect[rect] goodresp(nkeep) goodfevd(nkeep)
declare vector ik a(nvar) ones(nvar)
declare series[rect] irfsquared
compute ones=%const(1.0)
*
compute accept=0
infobox(action=define,progress,lower=1,upper=n1) $
   "Monte Carlo Integration"
do draws=1,n1
   *
   * Make a draw from the posterior for the VAR and compute its impulse
   * responses.
   *
   compute sigmad  =%ranwisharti(svt,wishdof)
   compute swish   =%decomp(sigmad)
   compute betau   =%ranmvkron(swish,sxx)
   compute betadraw=betaols+betau
   compute %modelsetcoeffs(varmodel,betadraw)
   impulse(noprint,model=varmodel,decomp=swish,results=impulses,steps=nstep)
   gset irfsquared 1 1     = %xt(impulses,t).^2
   gset irfsquared 2 nstep = irfsquared{1}+%xt(impulses,t).^2
   *
   * Do the subdraws over the unit sphere. These give the weights on the
   * orthogonal components.
   *
   do subdraws=1,n2
      compute a=%ransphere(nvar)
      *
      * Check that the responses have the correct signs for steps 1 to
      * KMAX+1 (+1 because in the paper, the steps used 0-based
      * subscripts, rather than the 1 based subscripts used by RATS).
      *
      if UhligAccept(a,KMIN,KMAX+1,||+2,+4,-5||)==0
         goto reject
      *
      * This is an accepted draw. Copy the information out. If we have
      * enough good ones, drop out of the overall loop.
      *
      compute accept=accept+1
      dim goodresp(accept)(nstep,nvar) goodfevd(accept)(nstep,nvar)
      ewise goodresp(accept)(i,j)=ik=%xt(impulses,i)*a,ik(j)
      ewise goodfevd(accept)(i,j)=$
         ik=(irfsquared(i)*(a.^2))./(irfsquared(i)*ones),ik(j)
      if accept>=nkeep
         break
   :reject
   end do subdraws
   if accept>=nkeep
      break
   infobox(current=draws)
end do draws
infobox(action=remove)
*
* Post-processing. Graph the mean of the responses along with the 16%
* and 84%-iles.
 clear upper lower resp

dec vect[series] targets(nvar)
do i=1,nvar
set targets(i) 1 nstep = 0.0
smpl 1 accept
do k=1,nstep
set work = goodresp(t)(k,i)
compute frac=%fractiles(work,||.50||)
compute targets(i)(k)=frac(1)
end do k
end do i
*
dec vect g(nvar-1)
compute g=%fill(nvar-1,1,1.0)
nonlin g
dec real gap
find(trace) min gap
compute a=%stereo(g)
compute gap=0.0
do k=1,nstep
compute gap=gap+%normsqr(%xt(targets,k)-%xt(impulses,k)*a)
end do k
end find
*
set resp 1 nstep = 0.0


*
spgraph(vfields=3,hfields=2,$
  footer="Figure 6. Impulse Responses with Pure-Sign Approach")
do i=1,nvar
   set resp 1 nstep = ik=%xt(impulses,t)*a,ik(i)
graph(ticks,number=0,picture="##.##",header="Impulse Responses for "+vl(i))
# resp 1 nstep
end do i

*
spgraph(done)

** in seperate graph
do i=1,nvar
   smpl 1 accept
   do k=1, nstep
   set work = goodresp(t)(k,i)
   compute frac=%fractiles(work, ||.16,.84||)
   compute lower(k)=frac(1)
   compute upper(k)=frac(2)
   end do k

   set resp 1 nstep = ik=%xt(impulses,t)*a,ik(i)


   *
   smpl 1 nstep
    graph(ticks,number=0) 3
   # resp
# upper / 2
   # lower / 2
end do i