When I run the following code I get the following error message. I am using rats version 9. I also attached the data for your check.
## MAT14. Non-invertible Matrix. Using Generalized Inverse for SYMMETRIC.
Code: Select all
************
************ This program, coupled with the procedure VARTVPKSC.src estimates the
************ VAR with time-varying parameters and stochastic volatility of
************ Giorgio Primiceri, (2005), "Time Varying Structural Vector Autoregressions and Monetary Policy,"
************ Review of Economic Studies, 72, 821-852.
************ The program and procedure were written by Todd Clark, todd.e.clark@kc.frb.org
dis %ratsversion()
dis %dateandtime()
comp dropexpl = 0 ;* 0 to not drop explosive draws in MCMC
comp ndraws = 8000 ;* number of draws retained. total draws = burndraws + ndraws*skipinterval
comp burndraws = 5000
comp skipinterval = 2 ;* set this to 1 to not skip any draws in what is retained
comp usemedians = 1 ;* 1 to use posterior medians in most charts, 0 to use posterior means
comp commonscale = 1 ;* 1 to impose common scale in impulse response plots, 0 to not
grparm header 18 subheader 18 keylabeling 16
env nowshowgraphs gsave="repPrimiceri_*.rgf" gformat=rgf
comp nvar = 6 ;* number of variables in VAR
comp fixlags = 1 ;* fixed lag order in model
comp styr = 1985
cal(m) styr:2
comp stsmpl = styr:2
comp endsmpl = 2015:4 ;* "endsmpl" is the last observation period.
all stsmpl endsmpl
smpl stsmpl endsmpl
comp seedval = 44*48
seed seedval
********************************
******************************** reading in Koop's data
********************************
open data 2015_expass.xlsx
data(for=xlsx,org=col) / loil lindpro ler lm lwpi lcpi
close
** for storing variable names, allocate label vector and order position of real variable
dec vec[string] varlabels(nvar)
dec vec[ser] y(nvar)
comp k = 0
dofor i = loil lindpro ler lm lwpi lcpi
comp k = k+1
set y(k) = i{0}
comp varlabels(k) = %l(i)
end do i
********************************
******************************** setting data sample, and weighting energy inflation
********************************
comp prenobs = 72
comp stpt = 1991:2
comp endpt = endsmpl
********************************
******************************** setting dimensions and prior stuff
******************************** Priors are same as in Primiceri.
sou(noecho) VARTVPKSC.src
** prior on VAR coefs (not used in this example)
dec vec meanfirstlag(nvar)
comp meanfirstlag = %const(0.)
** prior on time variation in coefficients (mean and IW df)
comp Q_priordf = prenobs
comp Q_priorscale = .0001
** prior on time variation in each equation's A coefficients (mean and IW df)
** (equation 1 stuff here not actually used)
dec vec[int] S_priordf(nvar)
do i = 1,nvar
comp S_priordf(i) = i
end do i
comp S_priorscale = .01
** prior on time variation in log volatilities, which are log st devs (mean and IW df)
comp priordfPhi = nvar+1
dec symm muPhi(nvar,nvar)
comp muPhi = %mscalar(0.0001)
** prior variance of initial value of volatility (mean set by procedure)
dec symm OmegaDelta0(nvar,nvar)
comp OmegaDelta0 = %mscalar(1.)
********************************
******************************** OLS estimation, full and pre-sample, just as checks
********************************
system(model=varmodel)
variables y
lags 1 to fixlags
det constant
end(system)
************ full sample
estimate(noftests) stpt endpt
comp compcoef = %modelcompanion(varmodel)
eigen(cvalues=cxvals,sort=absval) compcoef
dis 'largest eigen value in full sample = ' %real(%cabs(cxvals(1)))
comp OLSsigma = %sigma ;* will use this below to set fixed, common size of shocks
************ training sample
estimate(noftests) stpt-prenobs stpt-1
comp compcoef = %modelcompanion(varmodel)
eigen(cvalues=cxvals,sort=absval) compcoef
dis 'largest eigen value in training sample = ' %real(%cabs(cxvals(1)))
********************************
******************************** estimation, along with reporting of a limited set of posterior statistics
********************************
comp meanfirstlag = %const(.0)
comp cbar = .001
@VARTVPKSC(prmean,prinit,inclconst,doRsq,nodropexpl,noinformprior) y stpt endpt prenobs fixlags ndraws burndraws skipinterval meanfirstlag Q_priordf Q_priorscale $
S_priordf S_priorscale muDelta0 OmegaDelta0 muPhi priordfPhi cbar PiRes ARes Qres PhiRes LambdaRes stdevRes SRes RsqRes
dis %dateandtime()
********************************
******************************** extraction of various results of interest for intercepts and VAR slope coefs
********************************
********** computing means of each variable implied by model (all draws)
comp ncoef = fixlags*nvar+1 ;* # of coefficients in each equation
dec vec[rec] means_all(nvar)
do i = 1,nvar
dim means_all(i)(ndraws,endpt)
end do i
dec vec pivec(nvar*ncoef) compmeans compinterc
dec rec compcoef idenmatrix
** just to set some dimensions and avoid repetitive calcs
comp compcoef = %modelcompanion(varmodel)
comp idenmatrix = %identity(%rows(compcoef))
comp compinterc = %zeros(%rows(compcoef),1)
do draws = 1,ndraws
do vtime = stpt,endpt
** assigning model coefs, and set up companion form coefs
do i = 1,nvar*ncoef
comp pivec(i) = pires(draws,i)(vtime)
end do i
comp %modelsetcoeffs(varmodel,%vectorect(pivec,ncoef))
comp compcoef = %modelcompanion(varmodel)
** compute means
do i = 1,nvar
comp compinterc(i) = pivec(i*ncoef)
end do i
comp compmeans = inv(idenmatrix - compcoef)*compinterc
** storage of computed means
do i = 1,nvar
comp means_all(i)(draws,vtime) = compmeans(i)
end do i
end do vtime
end do draws
********** forming sums of AR coefs (all draws)
dec rec[rec] sar(nvar,nvar) ;* rows: equations; columns: variables
do i = 1,nvar
do j = 1,nvar
dim sar(i,j)(ndraws,endpt)
end do j
end do i
comp xpi = %zeros(fixlags,1) ;* for pulling coefs for variable i in equation j
comp [vec] pikeep = %zeros(ncoef*nvar,1)
do i = 1,ndraws
do j = stpt,endpt
do ii = 1,ncoef*nvar
comp pikeep(ii) = pires(i,ii)(j)
end do ii
do ii = 1,nvar
do jj = 1,nvar
comp xpi = %xsubvec(pikeep,(ii-1)*ncoef+(jj-1)*fixlags+1,(ii-1)*ncoef+jj*fixlags)
comp sar(ii,jj)(i,j) = %sum(xpi)
end do jj
end do ii
end do j
end do i
********** forming posterior means and 70% intervals for sums of AR coefs
smpl stpt endpt
dec rec[ser] sarmean(nvar,nvar) sar15(nvar,nvar) sar85(nvar,nvar)
clear sarmean sar15 sar85
smpl 1 ndraws
do vtime = stpt,endpt
do i = 1,nvar
do j = 1,nvar
set(scratch) statser = sar(i,j)(t,vtime)
stats(noprint,fractiles) statser
comp [vec] frac = %fractiles(statser,||.15,.5,.85||)
comp sarmean(i,j)(vtime) = %if(usemedians==1,%median,%mean)
comp sar15(i,j)(vtime) = frac(1)
comp sar85(i,j)(vtime) = frac(3)
end do j
end do i
end do vtime
********** forming posterior means and 70% interval for intercepts of each equation, and implied means
smpl stpt endpt
dec vec[ser] intercept_mean(nvar) intercept_15(nvar) intercept_85(nvar)
dec vec[ser] means_mean(nvar) means_15(nvar) means_85(nvar)
clear intercept_mean intercept_15 intercept_85 means_mean means_15 means_85
smpl 1 ndraws
do vtime = stpt,endpt
do i = 1,nvar
*
set(scratch) statser = PiRes(t,i*ncoef)(vtime)
stats(noprint,fractiles) statser
comp [vec] frac = %fractiles(statser,||.15,.5,.85||)
comp intercept_mean(i)(vtime) = %if(usemedians==1,%median,%mean)
comp intercept_15(i)(vtime) = frac(1)
comp intercept_85(i)(vtime) = frac(3)
*
set(scratch) statser = means_all(i)(t,vtime)
stats(noprint,fractiles) statser
comp [vec] frac = %fractiles(statser,||.15,.5,.85||)
comp means_mean(i)(vtime) = %if(usemedians==1,%median,%mean)
comp means_15(i)(vtime) = frac(1)
comp means_85(i)(vtime) = frac(3)
*
end do i
end do vtime
********************************
******************************** extraction of various results of interest relating to stochastic volatility
********************************
********** forming posterior means and 70% intervals for elements of A
smpl stpt endpt
dec vec[ser] Amean A15 A85
if nvar>1
{
dim Amean(nvar*(nvar-1)/2) A15(nvar*(nvar-1)/2) A85(nvar*(nvar-1)/2)
clear Amean A15 A85
}
endif
if nvar>1
{
smpl 1 ndraws
do vtime = stpt,endpt
do i = 1,nvar*(nvar-1)/2
set(scratch) statser = ARes(t,i)(vtime)
stats(noprint,fractiles) statser
comp [vec] frac = %fractiles(statser,||.15,.5,.85||)
comp Amean(i)(vtime) = %if(usemedians==1,%median,%mean)
comp A15(i)(vtime) = frac(1)
comp A85(i)(vtime) = frac(3)
end do i
end do vtime
}
endif
********** forming posterior means and 70% intervals for structural and reduced form shock st devs
smpl stpt endpt
dec vec[ser] Deltamean Delta15 Delta85 stdevmean stdev15 stdev85
dim Deltamean(nvar) Delta15(nvar) Delta85(nvar) stdevmean(nvar) stdev15(nvar) stdev85(nvar)
clear Deltamean Delta15 Delta85 stdevmean stdev15 stdev85
smpl 1 ndraws
do vtime = stpt,endpt
do i = 1,nvar
set(scratch) statser = sqrt(LambdaRes(t,i)(vtime))
stats(noprint,fractiles) statser
comp [vec] frac = %fractiles(statser,||.15,.5,.85||)
comp Deltamean(i)(vtime) = %if(usemedians==1,%median,%mean)
comp Delta15(i)(vtime) = frac(1)
comp Delta85(i)(vtime) = frac(3)
*
set(scratch) statser = stdevRes(t,i)(vtime)
stats(noprint,fractiles) statser
comp [vec] frac = %fractiles(statser,||.15,.5,.85||)
comp stdevmean(i)(vtime) = %mean ;* just force mean, to because that is what Koop reports, even though he reports medians of imp resp
*comp stdevmean(i)(vtime) = %if(usemedians==1,%median,%mean)
comp stdev15(i)(vtime) = frac(1)
comp stdev85(i)(vtime) = frac(3)
end do i
end do vtime
********** forming posterior means and 70% intervals for R2
smpl stpt endpt
dec vec[ser] r2mean(nvar) r215(nvar) r285(nvar)
clear r2mean r215 r285
smpl 1 ndraws
do vtime = stpt,endpt
do i = 1,nvar
set(scratch) statser = RsqRes(t,i)(vtime)
stats(noprint,fractiles) statser
comp [vec] frac = %fractiles(statser,||.15,.5,.85||)
comp r2mean(i)(vtime) = %if(usemedians==1,%median,%mean)
comp r215(i)(vtime) = frac(1)
comp r285(i)(vtime) = frac(3)
end do i
end do vtime
********************************
******************************** dimensioning stuff used below in impulse responses, etc.
********************************
comp nstep = 16 ;* number of impulse response steps
dec vec shocksize(nvar) scalefact(nvar)
comp shocksize = %sqrt(%xdiag(OLSsigma)) ;* fixed shock size considered in IR's, based on full sample shock sizes
comp [vec[int]] impdates = ||1994:1,2001:4,2015:1||
comp nimp = %rows(impdates) ;* number of periods for which I compare impulse responses
dec vec[str] datestring(nimp)
do i = 1,nimp
comp datestring(i) = %datelabel(impdates(i))
end do i
dec rec[rec] responses(nimp,ndraws) ;* matrix for storing draws of impulse responses based on recursive
comp [symm] Lambdamat = %zeros(nvar,nvar)
dec symm Sigma
comp Amat = %identity(nvar)
dec vec pivec(nvar*ncoef)
********************************
******************************** computing impulse responses, using draws of VAR coefs and A (fixed shock size)
********************************
do count = 1,nimp
comp vtime = impdates(count)
do draws = 1,ndraws
dim responses(count,draws)(nvar*nvar,nstep)
** setting A and Sigma matrices
comp pos = 0
do i = 2,nvar
do j = 1,i-1
comp pos = pos + 1
comp Amat(i,j) = ARes(draws,pos)(vtime)
end do j
end do i
comp invAmat = inv(Amat)
do i = 1,nvar
comp Lambdamat(i,i) = LambdaRes(draws,i)(vtime)
end do i
comp Sigma = %mqform(Lambdamat,tr(invAmat))
comp swish = %decomp(Sigma)
** assigning model coefs, and computing impulse responses
do i = 1,nvar*ncoef
comp pivec(i) = pires(draws,i)(vtime)
end do i
comp %modelsetcoeffs(varmodel,%vectorect(pivec,ncoef))
impulse(model=varmodel,decomp=swish,steps=nstep,results=impulses,noprint)
** storage in "impulses": series i,j is the response of variable i to a shock to j
** storing CHOLESKI-BASED responses. RESPONSES, cols 1-nvar: response of each variable to shock 1; cols 2-nvar: resp of each var to shock 2; etc
ewise responses(count,draws)(i,j)=impulses(%clock(i,nvar),(i-1)/nvar+1)(j)
end do draws
************ now normalize the responses for this period to make the average shock sizes the same across time (= full sample OLS size)
** compute scaling factors from median shock sizes, Choleski
do i = 1,nvar
set statser 1 ndraws = responses(count,t)((i-1)*nvar+i,1)
stats(fractiles,noprint) statser 1 ndraws
comp scalefact(i) = shocksize(i)/%median
end do i
** now do the rescaling
do draws = 1,ndraws
do k = 1,nstep
do i = 1,nvar
do j = 1,nvar
comp responses(count,draws)((i-1)*nvar+j,k) = responses(count,draws)((i-1)*nvar+j,k)*scalefact(i)
end do j
end do i
end do k
end do draws
end do count
********************************
******************************** charts
********************************
********** global settings for charts
smpl stpt endpt
if usemedians==1
{
comp [vec[string]] mykey_full = ||'median','15%ile','85%ile'||
comp [vec[string]] mykey_short = ||'median'||
}
else
{
comp [vec[string]] mykey_full = ||'mean','15%ile','85%ile'||
comp [vec[string]] mykey_short = ||'mean'||
}
endif
if nvar==5.or.nvar==6
{
comp vf = 2
comp hf = 3
}
else if nvar==3.or.nvar==4
{
comp vf = 2
comp hf = 2
}
endif
if nvar==7.or.nvar==6
{
comp vf_a = 2
comp hf_a = 3
}
else if nvar==5.or.nvar==4
{
comp vf_a = 2
comp hf_a = 2
}
else
{
comp vf_a = 2
comp hf_a = 1
}
endif
********** charts: sums of coefs -- means and percentiles
do i = 1,nvar
comp subheader = 'equation for ' + varlabels(i)
spgraph(vfields=vf,hfields=hf,header='TVP-VAR: sums of coefficients',subheader=subheader)
do j = 1,nvar
comp header = 'lags of variable: ' + varlabels(j)
graph(dates,header=header,key=below,klab=mykey_full) 3
# sarmean(i,j) / 1
# sar15(i,j) / 2
# sar85(i,j) / 2
end do j
spgraph(done)
end do i
********** charts: intercepts -- means and percentiles, and then just equation means
list meancard = 1 to nvar
spgraph(vfields=vf,hfields=hf,header='TVP-VAR: intercepts for each equation')
do i = 1,nvar
comp header = 'equation for ' + varlabels(i)
graph(dates,header=header,key=below,klab=mykey_full) 3
# intercept_mean(i) / 1
# intercept_15(i) / 2
# intercept_85(i) / 2
end do i
spgraph(done)
spgraph(vfields=vf,hfields=hf,header='TVP-VAR: means or medians for each variable')
do i = 1,nvar
comp header = varlabels(i)
graph(dates,header=header,key=below,klab=mykey_full) 3
# means_mean(i) / 1
# means_15(i) / 2
# means_85(i) / 2
end do i
spgraph(done)
graph(dates,header='TVP-VAR: means or medians for all variables together',key=below) nvar
cards means_mean(meancard)
********** charts: A coefs and percentiles
if nvar>1
{
comp pos = 0
do i = 2,nvar
comp header = 'TVP-VAR: A coefficients, equation for ' + varlabels(i)
spgraph(vfields=vf_a,hfields=hf_a,header=header)
do j = 1,i-1
comp pos = pos + 1
comp header = 'coef on ' + varlabels(j)
graph(dates,header=header,key=below,klab=mykey_full) 3
# Amean(pos) / 1
# A15(pos) / 2
# A85(pos) / 2
end do j
spgraph(done)
end do i
}
endif
********** charts: residual standard deviations, structural vs reduced form
spgraph(vfields=vf,hfields=hf,header='TVP-VAR: residual standard deviations')
do i = 1,nvar
set reducedform = stdevmean(i){0}
set structural = Deltamean(i){0}
comp header = 'equation for ' + varlabels(i)
graph(dates,header=header,key=below,klab=||'reduced form','structural'||) 2
# reducedform
# structural
dis '****************** reduced form and structural residual standard deviations of ' varlabels(i)
print / reducedform structural
end do i
spgraph(done)
********** charts: residual standard deviations, reduced form with 70 pct bands
spgraph(vfields=vf,hfields=hf,header='TVP-VAR: reduced form st devs and 70pct interval')
do i = 1,nvar
comp header = 'equation for ' + varlabels(i)
graph(dates,header=header,key=below,klab=mykey_full) 3
# stdevmean(i)
# stdev15(i) / 2
# stdev85(i) / 2
end do i
spgraph(done)
********** impulse responses based on Choleski ID: displaying estimates and creating one chart comparing all inflation responses over time
dec vec[series] upper(nvar) lower(nvar) resp(nvar)
smpl 1 nstep
clear upper lower resp
*** change some graph settings to make things more readable in the impulse response charts (therefore, we want response charts to come last)
grparm(bold) hlabel 16 matrixlabels 14 header 18
grparm axislabel 36
do count = 1,nimp
comp vtime = impdates(count)
dis ''
dis '************************************************'
dis '************************************************ impulse responses for' %datelabel(vtime)
dis '************************************************'
comp header = 'Impulse responses: ' + %datelabel(vtime)
do i = 1,nvar
comp subheader = 'shock to ' + varlabels(i)
do j=1,nvar
comp [string] chartheader = 'response of ' + varlabels(j)
smpl 1 ndraws
do k=1,nstep
set(scratch) work 1 ndraws = responses(count,t)((i-1)*nvar+j,k)
stats(noprint,fractiles) work 1 ndraws
compute frac=%fractiles(work,||.15,.85||)
compute lower(j)(k)=frac(1)
compute upper(j)(k)=frac(2)
compute resp(j)(k)=%if(usemedians==1,%median,%mean)
end do k
end do j
spgraph(done)
dis '*********** posterior means/medians of responses to shock in ' varlabels(i)
print(nodates) 1 nstep resp
dis '*********** 15th percentile of responses to shock in ' varlabels(i)
print(nodates) 1 nstep lower
dis '*********** 85th percentile of responses to shock in ' varlabels(i)
print(nodates) 1 nstep upper
end do i
end do count
** now do plot comparing responses over time to inflation shock
smpl 1 nstep
if usemedians==1
comp subheader = '(posterior medians)'
else
comp subheader = '(posterior means)'
endif
********** charts: impulse responses, all for a given time period on a single page
do count = 1,nimp
comp vtime = impdates(count)
comp header = 'Impulse responses: ' + %datelabel(vtime)
spgraph(header=header,subheader='(point estimates with 70% bands)',xpos=both,xlab=varlabels, $
ylab=varlabels,vlab="Responses of",hlab="Shocks to",vfields=nvar,hfields=nvar)
*
* Because we want a common scale for all responses of a single variable, we need
* to do all the calculations for a full row of graphs first.
*
* dec vect[series] upper(nvar) lower(nvar) resp(nvar)
do i=1,nvar
compute minlower=maxupper=0.0
smpl 1 ndraws
do j=1,nvar
clear lower(j) upper(j) resp(j)
do k=1,nstep
set(scratch) work 1 ndraws = responses(count,t)((j-1)*nvar+i,k)
stats(noprint,fractiles) work 1 ndraws
compute frac=%fractiles(work,||.15,.85||)
compute lower(j)(k)=frac(1)
compute upper(j)(k)=frac(2)
compute resp(j)(k)=%if(usemedians==1,%median,%mean)
end do k
compute maxupper=%max(maxupper,%maxvalue(upper(j)))
compute minlower=%min(minlower,%minvalue(lower(j)))
end do j
*
smpl 1 nstep
do j=1,nvar
if commonscale==1
{
graph(ticks,min=minlower,max=maxupper,number=0) 3 j i
# resp(j)
# upper(j) / 2
# lower(j) / 2
}
else
{
graph(ticks,number=0) 3 j i
# resp(j)
# upper(j) / 2
# lower(j) / 2
}
end do j
end do i
*
*
spgraph(done)
end do count
dis %dateandtime()