******************
*Application Bojrlond 2009
******************
*SVAR monetary
***

open data data.rat
calendar(q) 1994
data(format=rats) 1994:1 2012:4 REER DPG TXINTER PIB IPC
*
* Data are already in logs. We need to scale by 100.

****************************************
****transformation
****************************************

*set c  = log(c)
*set dc  = c-c(-1)
*set in = in*100.0
*
set lreer = log(reer)
set ldpg = log(dpg)
set txinter = txinter
set lpib = log(pib)
set lipc = log(ipc)
*
diff lipc / dlipc
esmooth(trend=select,seasonal=select,smooth=sldpg,$
initial=start) ldpg
**gaph depenses desaisonalisees
graph(footer="Figure 2",$
  key=below,klabels=||"Depenses publiques globales"||) 1
#sldpg
**

FILTER(TYPE=HP,TUNING=1600) LPIB 1994:01 2012:04 LPIBt
set lpibhp = lpib - lpibt

FILTER(TYPE=HP,TUNING=1600) txinter 1994:01 2012:04 txintert
set txinterhp = txinter - txintert

FILTER(TYPE=HP,TUNING=1600) sldpg 1994:01 2012:04 sldpgt
set SLDPGhp = sldpg - sldpgt

FILTER(TYPE=HP,TUNING=1600) Lreer 1994:01 2012:04 Lreert
set LREERHP  = Lreer - Lreert

***********
*Estimation
***********

compute nsteps = 30
*
system(model=Ben)
var dlipc lpibhp txinterhp SLDPGhp LREERHP
lags 1 to 4
det constant
end(system)
estimate(resids=resids) 1998:1 2012:4
compute vsigma=%sigma


********************************************
* Model with short and long-run restrictions
********************************************

*
dec rect lr(5,5) sr(5,5)

input lr
. 0 . . .
. . 0 . .
0 . . . .
. . 0 . 0
. . . . .
input sr
. 0 . . 0
. . 0 . .
. 0 . . .
. . 0 . .
. . . . .
dec vect[strings] shocklabels varlabels
compute shocklabels=||"Inflation","supply","Monetary policy","budgetary policy","foreign"||
compute varlabels=||"dlipc","dlpib","dtxinter","dsldpg","dlreer"||
*
@ShortAndLong(sr=sr,lr=lr,masum=inv(%varlagsums)) %sigma f
compute fd=%xdiag(f)
ewise f(i,j)=f(i,j)/fd(j)
*
* Point estimates
*
@varirf(model=ben,steps=nsteps,factor=f,page=byshock,$
  shocks=shocklabels,varlabels=varlabels,accumulate=||1,2,3,4,5||)


***************************
* Monte Carlo integration
***************************
*
procedure SRLRDoDraws
*
option model     model
option integer   draws   10000
option integer   steps   40
option vect[int] accum   ||1||
option rect      lr
option rect      sr
*
local integer nvar
local rect    fxx fwish fsigmad factor
local integer wishdof
local rect    betaols betau betadraw
local vect    ix
local symm    sigmad
local integer i j draw
*
if .not.%defined(model) {
   disp "###SRLRDoDraws(MODEL=model name,other options)"
   return
}
compute nvar=%modelsize(model)
*
* Standard setup for drawing from an OLS VAR
*
compute fxx    =%decomp(%xx)
compute fwish  =%decomp(inv(%nobs*%sigma))
compute wishdof=%nobs-%nreg
compute betaols=%modelgetcoeffs(model)
*
local rect[series] impulses(nvar,nvar)
*
* These are global variables
*
declare vect[rect]   %%responses(draws)
*
infobox(action=define,progress,lower=1,upper=draws) "Monte Carlo Integration"
do draw=1,draws
   if %clock(draw,2)==1 {
      compute sigmad  =%ranwisharti(fwish,wishdof)
      compute fsigmad =%decomp(sigmad)
      compute betau   =%ranmvkron(fsigmad,fxx)
      compute betadraw=betaols+betau
   }
   else
      compute betadraw=betaols-betau
   compute %modelsetcoeffs(model,betadraw)
   *
   * Compute the short-and-long-run factor using the recalculated lag
   * sums of the model.
   *
   @ShortAndLong(lr=lr,sr=sr,masum=inv(%modellagsums(model))) sigmad factor
   *
   * The rescaling done in the next two lines is specific to this paper
   *
   compute fd=%xdiag(factor)
   ewise factor(i,j)=factor(i,j)/fd(j)
   *
   impulse(noprint,model=model,factor=factor,results=impulses,steps=steps)
   *
   * Accumulate the responses as requested
   *
   do i=1,%rows(accum)
      do j=1,nvar
         acc impulses(accum(i),j) 1 steps
      end do j
   end do i
   *
   * Store the impulse responses
   *
   dim %%responses(draw)(nvar*nvar,steps)
   ewise %%responses(draw)(i,j)=ix=%vec(%xt(impulses,j)),ix(i)
   infobox(current=draw)
end do draw
infobox(action=remove)
*
* Restore the original coefficients
*
compute %modelsetcoeffs(model,betaols)
*
end SRLRDoDraws
*
@SRLRDoDraws(accum=||2,3,4||,steps=nsteps,model=ben,lr=lr,sr=sr)
@MCProcessIRF(model=ben,percentiles=||.16,.84||,center=median,lower=lower,upper=upper,irf=irf)
*
spgraph(vfields=2,footer=$
  "Figure 2 Impulse Responses to a Exchange rate shock and expenditure Shock")
  do block=1,2
     if block==1
        compute shockcolumn=5,shockheader="Exchange rate shock"
     else
        compute shockcolumn=4,shockheader="expenditure Shock"


spgraph(vfields=5,hfields=5,header=shockheader)
     graph(row=1,col=1,header=varlabels(5),nodates) 3
     # irf(5,shockcolumn)   / 1
     # lower(5,shockcolumn) / 2
     # upper(5,shockcolumn) / 2

graph(row=1,col=2,header=varlabels(4),nodates) 3
     # irf(4,shockcolumn)   / 1
     # lower(4,shockcolumn) / 2
     # upper(4,shockcolumn) / 2

graph(row=4,col=1,header=varlabels(3),nodates) 3
     # irf(3,shockcolumn)   / 1
     # lower(3,shockcolumn) / 2
     # upper(3,shockcolumn) / 2


     graph(row=2,col=1,header=varlabels(2),nodates) 3
     # irf(2,shockcolumn)   / 1
     # lower(2,shockcolumn) / 2
     # upper(2,shockcolumn) / 2

     graph(row=2,col=2,header=varlabels(1),nodates) 3
     # irf(1,shockcolumn)   / 1
     # lower(1,shockcolumn) / 2
     # upper(1,shockcolumn) / 2
spgraph(done)
  end do block
spgraph(done)