clear(all)
environnement noecho
calendar(panelobs=42,a) 1970:1
allocate 20//42
open data database1.xls
data(org=obs,format=xls)


compute nbvar=3;   * number of variables in the PVAR model
*=====================================================================
* List of endogenous variables in the PVAR model
dec vect[strings]  ylabel(nbvar)  varslabel
dec vector[int]    endo(nbvar)    varendo

compute endo(1)=X1;       compute ylabel(1)= %l(endo(1));
compute endo(2)=X2;       compute ylabel(2)= %l(endo(2));
compute endo(3)=X3;       compute ylabel(3)= %l(endo(3));

compute varendo = %rlempty()
do i=1,nbvar
   compute varendo=%rladdone(varendo,endo(i))
end do i
compute nperiods = %panelobs();
sta(noprint) country; compute nobscountry = %nobs;
*=====================================================================
* List of parameters in the PVAR model
compute nblags=2;                      * number of lags in the PVAR model
compute nby = nperiods;                * number of time periods per country
compute nbc = nobscountry/nperiods;    * number of countries
compute nstep=11;          * number of forecats steps
compute ndraws=1000;       * number of replications in Monte Carlo simulations
compute f1=0.10;           * lower percentile of the impluse-response distribution in the figures
compute f2=1-f1;           * upper percentile of the impluse-response distribution in the figures
compute gg=1;              * position of the exogenous variable
*=====================================================================
* Country fixed effect, country-specific linear and quadractic trends and time fixed effect
set ltrend = %period(t)
dec vector[series] dummies(nbc) trend(nbc) trend2(nbc) time(nby)
do i=1,nbc
	set dummies(i) = %indiv(t)==i;                    * country fixed effect
	set trend(i)   = (%indiv(t)==i)*ltrend;           * country-specific linear trend
     set trend2(i)  = (%indiv(t)==i)*ltrend*ltrend;    * country-specific quadratic trend
end do i
do i=1,nby
set time(i)  = (%period(t)==i);                        * time fixed effect
end do i

smpl

*=====================================================================
* VAR system
system(model=pvarmodel1)
variables varendo
lags 1 to nblags
kfset XXX
det  trend dummies
end(system)
estimate(noprint,outsigma=vmat)





*=====================================================================
* Montecarlo Simulation for IRFs bands
declare rect sxx swish betaols betadraw
declare symm sigmad
declare vect[rect] responses(ndraws)
declare rect[series] impulses(nbvar,nbvar)
compute sxx    =%decomp(XXX);                compute svtr   =tr(%decomp(vmat));
compute betaols=%modelgetcoeffs(pvarmodel1); compute ncoef  =%rows(sxx);
compute wishdof=%nobs-ncoef;
declare rect ranc(ncoef,nbvar)
infobox(action=define,progress,lower=1,upper=ndraws) 'Monte Carlo Integration'
do draws = 1,ndraws
   if %clock(draws,2)==1 {
      compute sigmad  =%mqform(%nobs*inv(%ranwishart(nbvar,wishdof)),svtr)
      compute swish   =%decomp(sigmad)
	 compute ranc    =%ran(1.0)
      compute betau   =sxx*ranc*tr(swish)
      compute betadraw=betaols+betau
   }
   else
      compute betadraw=betaols-betau
   compute %modelsetcoeffs(pvarmodel1,betadraw)
   impulse(noprint,model=pvarmodel1,decomp=swish,results=impulses) nbvar nstep
   dim responses(draws)(nbvar*nbvar,nstep)
   ewise responses(draws)(i,j)=impulses((i-1)/nbvar+1,%clock(i,nbvar))(j)
   infobox(current=draws)
end do draws
infobox(action=remove)

* FIGURES IRFs (shock in fiscal)
*=====================================================================
dec vect[series] upper(nbvar) lower(nbvar) resp(nbvar)
dec vect[series] irfm(nbvar) irfu(nbvar) irfl(nbvar) cirfm(nbvar)

do j=1,nbvar
  compute varslabel = varslabel+ylabel(j)
end do


compute shortheadervar = varslabel(1)+' '+varslabel(2)+' '+varslabel(3)
compute headervar = 'IRF to X1 Shock ('+ shortheadervar +')'


grparm(italics)
spgraph(xpos=upper, ylab=ylabel,vfields=nbvar,hfields=1,header=headervar)
smpl 1 ndraws
clear resp(gg)
do k=1,nstep
set work 1 ndraws = responses(t)((gg-1)*nbvar+gg,k)
compute frac=%fractiles(work,||f1,f2||)
compute resp(gg)(k)=%avg(work)
end do k
compute ggg=resp(gg)(1)
do i=1,nbvar
   compute minlower=maxupper=0.0
   smpl 1 ndraws
   clear lower(gg) upper(gg) resp(gg)
      do k=1,nstep
         set work 1 ndraws = responses(t)((i-1)*nbvar+gg,k)
         compute frac=%fractiles(work,||f1,f2||)
         compute lower(gg)(k)=frac(1)/ggg
         compute upper(gg)(k)=frac(2)/ggg
         compute resp(gg)(k)=%avg(work)/ggg
     end do k
	compute maxupper=%max(maxupper,%maxvalue(upper(gg)))
	compute minlower=%min(minlower,%minvalue(lower(gg)))
      	smpl 1 nstep
set irfm(i) = resp(gg); set irfu(i) = upper(gg); set irfl(i) = lower(gg);
graph(min=minlower,max=maxupper,number=0,samesize) 3 1 i
      	# resp(gg)
      	# upper(gg) / 2
      	# lower(gg) / 2
end do i


spgraph(done)

do j = 1,nbvar
accumulate irfm(j) / cirfm(j)
end do j