Help with PVAR code
Posted: Tue Sep 12, 2017 9:39 am
I have adapted this code to run with my data (sorry the code is not fully cleaned). I am trying to get the impulse response graphs to appear one at the time. Also the code seems to print black and white and I would like them in color, please. Any help please?
Code: Select all
************* ESTIMATION OF PVAR FOR 18 COUNTRIES OVER 1970-2008 PERIOD ***********
environment noecho
calendar(panelobs=39) 1970
allocate 18//2008:1
open data D:\temp\Rats\PVAR\crisis.xls
data(format=xls, org=columns)
*
*
COMPUTE LAGS=2 ;*Set number of logs
COMPUTE NVAR=5
COMPUTE NSTEP=11
COMPUTE NDRAWS=2000
COMPUTE N = 18 ;*Number of countries
compute YYY = 37 ;*Number of observations for each country
compute f1=0.05
compute f2=0.95
compute gg=1 ;*Position of Government Spending G
compute tt=2 ;*Position of Government Net Taxes NT
compute yy=3 ;*Position of Output Y
compute ratiogy=0.241 ;*Average Ratio G/Y
compute rationty=0.217 ;*Average Ratio NT/Y
*
*
dec vect[strings] ylabel(nvar+2)
compute ylabel(gg)= 'Exp'
compute ylabel(tt)= 'Inv'
compute ylabel(yy)= 'Crisis'
compute ylabel(4)= 'Tax '
compute ylabel(5)= 'Debt'
compute ylabel(6)= 'Budget'
*
* this created country specific dummies and linear trend
*
SET LTREND = %PERIOD(T)
dec vector[series] dummies(N) time(YYY) trend(N) qtrend(N)
do i=1,N
set dummies(i) = %indiv(t)==i
set trend(i) = (%indiv(t)==i)*ltrend
set qtrend(i) = (%indiv(t)==i)*ltrend*ltrend
end do i
*
* time dummy
*
do i=1,YYY
set time(i) = (%period(t)==i)
end do i
*print / time
*
* ggflq{1}ggflq{2}
*
*********************************************************************************
* !!!!!!!!! change specification of the VAR !!!!!!!!!!!!!!!!!!!!!
**************************************************************************
system(model=varmodel)
variables g ntc y irl reer
lags 1 to lags
kfset xxx
det dummies time
end(system)
estimate(noprint,outsigma=vmat,RESIDUALS=RESVAR)
*set gshock = resvar(1)
*
* print / resvar
* The nonlinear parameters are put into a vector. This is not the
* most convenient way to code this for estimation (as it's harder
* to adjust the model), but it considerably cleans up the process
* of making draws.
*
* nfree is the number of free coefficients
*
**************************************************************************
* !!!!!!!!! change short-run identifying restrictions !!!!!!!!!!!!!!!!!
**************************************************************************
*
compute nfree=10
dec vect ax(nfree)
nonlin ax
*
dec frml[rect] afrml
frml afrml = || 1.0 ,0.0 ,0.0 ,0.0 ,0.0 |$
-ax(1) ,1.0 ,0.0 ,0.0 ,0.0 |$
-ax(10) ,-ax(3) ,1.0 ,0.0 ,0.0 |$
-ax(4) ,-ax(5) ,-ax(6) ,1.0 ,0.0 |$
-ax(7) ,-ax(8) ,-ax(9) ,-ax(2) ,1.0 ||
*
*
compute ax=%const(0.01)
*cvmodel(iters=20,method=simplex) vmat afrml
cvmodel(iters=300,method=bfgs,noprint) vmat afrml
*
declare rect sxx swish betaols betadraw
declare symm sigmad
*
declare rect sxx svtr swish betaols betadraw
declare symm sigmad
compute sxx =%decomp(xxx)
compute svtr =tr(%decomp(vmat))
compute betaols=%modelgetcoeffs(varmodel)
compute ncoef =%rows(sxx)
compute wishdof=%nobs-ncoef
*
declare rect ranc(ncoef,nvar)
*
declare vect[rect] responses(ndraws)
declare rect[series] impulses(nvar,nvar)
*
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(nvar,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(varmodel,betadraw)
impulse(noprint,model=varmodel,decomp=swish,results=impulses) nvar nstep
*
*
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)
*
**********************************************************************
* GRAPH
**********************************************************************
grparm(italics) axislabelling 25
grparm(nobold,italics) matrixlabelling 30
grparm header 30
*
spgraph(xpos=upper,vfields=3,hfields=2)
*
dec vect[series] upper(nvar) lower(nvar) resp(nvar)
*
*
smpl 1 ndraws
clear resp(gg)
do k=1,nstep
set work 1 ndraws = responses(t)((gg-1)*nvar+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,nvar
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)*nvar+gg,k)
*print(nodates) / work
compute frac=%fractiles(work,||f1,f2||)
compute lower(gg)(k)=frac(1)/(ggg*ratiogy)
compute upper(gg)(k)=frac(2)/(ggg*ratiogy)
compute resp(gg)(k)=%avg(work)/(ggg*ratiogy)
end do k
compute maxupper=%max(maxupper,%maxvalue(upper(gg)))
compute minlower=%min(minlower,%minvalue(lower(gg)))
smpl 1 nstep
display ylabel(i)
print / resp(gg) upper(gg) lower(gg)
graph(ticks,vticks=5,patterns,min=minlower,max=maxupper,number=0,samesize,picture='*.#',hlabel='years after shock',header=''+ylabel(i)) 3
# resp(gg) / 1
# upper(gg) / 6
# lower(gg) / 6
end do i
*
do i=6,6
compute minlower=maxupper=0.0
smpl 1 ndraws
clear lower(gg) upper(gg) resp(gg)
do k=1,nstep
set work 1 ndraws = -ratiogy*responses(t)((gg-1)*nvar+gg,k)+rationty*(responses(t)((tt-1)*nvar+gg,k)+2.1*responses(t)((yy-1)*nvar+gg,k))-(ratiogy-rationty)*responses(t)((yy-1)*nvar+gg,k)
compute frac=%fractiles(work,||f1,f2||)
compute lower(gg)(k)=frac(1)/(ggg*ratiogy)
compute upper(gg)(k)=frac(2)/(ggg*ratiogy)
compute resp(gg)(k)=%avg(work)/(ggg*ratiogy)
end do k
compute maxupper=%max(maxupper,%maxvalue(upper(gg)))
compute minlower=%min(minlower,%minvalue(lower(gg)))
smpl 1 nstep
display ylabel(i)
print / resp(gg) upper(gg) lower(gg)
graph(ticks,vticks=5,patterns,min=minlower,max=maxupper,number=0,samesize,picture='*.#',hlabel='years after shock',header=''+ylabel(i)) 3
# resp(gg) / 1
# upper(gg) / 6
# lower(gg) / 6
end do i
*
*
spgraph(done)
clear all
*