## MAT15. Subscripts Too Large or Non-Positive
The code is as follows:
Code: Select all
*
*
* This program computes average IRF for output, conditional on being in
* the upper or lower regime, averaging across initial conditions and then,
* for each initial setting, across bootstrapped residuals.
*
* Control parameters
*
comp nvar =5
comp horizon =15
*
* Number of bootstrap replications used in computing GIRF's
*
comp nkrep =500
*
* Change to upper=0 to get the lower regime
*
comp upper =1
*
comp [vector] shocksizes=||1.0,2.0||
comp [vector] shocksigns=||1.0,-1.0||
*
*
open data 20160516_ERPT.xls
calendar(q) 1980:1
data(format=xls,org=obs) 1980:01 2015:4
*
*
* These are the desired start and end of the sample to be used in the
* estimation. However, the transformations use the entire range.
*
compute sstart = 1982:1
compute send = 2015:4
*
*
* This is the exchange rate measure
*
set ner = dlner
*
dec vector[series] res(nvar) vres(nvar)
dec vector v(nvar) resmat(nvar)
dec vect[int] depvars(nvar)
dec rect[int] laglengths(nvar,nvar)
dec vect[equation] eqn(nvar)
dec vect[series] resup(nvar) resdn(nvar)
dec rect[frml] fitud(nvar,2)
dec vect[frml] tvarf(nvar) empty(nvar) tfrml(nvar) rfrml(nvar)
dec vect[series] bootu(nvar) data(nvar)
dec series[vect] bootuv bootres
*
dec vect[string] shortlabels(nvar) longlabels(nvar)
declare vector[labels] lab(nvar)
*
compute depvars =||dlpoilapsp,dlpricenpq,dlner,lgdpgap,dlcpi||
compute shortlabels=||"Oil","Rice","ER","YGAP","DLCPI"||
compute longlabels =||"World Oil Inflation","World Rice Inflation","Nominal Exchange Rate","Output Gap","Inflation Rate"||
compute lab =|| "Oil", "Rice", "ER", "YGAP","IR" ||
*
* set thresholds, mas, and delays
*
compute d = 1
*
* 15% window
*
dec vector macoeffs
compute malength = 2 ; comp thresh = 0.01582 ;* ner
*
* This generates the actual moving average of the data
*
filter(type=lagging,span=malength) dlner / nerthr
*
* This generates an equation (and from it a FRML) to compute the moving
* average of the data
*
dim macoeffs(malength)
comp macoeffs = %const(1./malength)
equation(identity,coeffs=macoeffs) threqn nerthr
# ner{0 to malength-1}
frml(equation=threqn,identity) thrfrml
*
* Lag lengths are allowed to differ among equations. These give the lag
* lengths with equation in a row and variable in a column.
*
input laglengths
4 4 4 4 4
4 4 4 4 4
4 4 4 4 4
4 4 4 4 4
4 4 4 4 4
*************************************************************************
*
* Dummy variables for the two regimes
*
set d1 = thrfrml{d}>thresh
set d2 = 1.-d1
*
compute maxlag=d+(malength-1)
inquire(reglist) rstart<<sstart rend<<send
# nerthr{d}
*
* Create an equation for each variable with the number of lags for each
* endogenous variable given by <<laglengths>>.
*
do i=1,nvar
compute [vect[integer]] rl=||constant||
do j=1,nvar
compute rl=%rladdlaglist(rl,depvars(j),%seq(1,laglengths(i,j)))
compute maxlag=%imax(maxlag,laglengths(i,j))
end do j
equation eqn(i) depvars(i)
# rl
end do i
*
compute rstart=sstart+maxlag,rend=send
do i=1,nvar
linreg(equation=eqn(i)) * rstart rend
frml empty(i) bootu(i) = 0.0
end do i
*
do i=1,nvar
disp
disp
disp shortlabels(i)+" regression in upper regime"
linreg(smpl=d1,equation=eqn(i),frml=fitud(i,1)) * rstart rend resup(i)
disp
disp
disp shortlabels(i)+" regression in lower regime"
linreg(smpl=d2,equation=eqn(i),frml=fitud(i,2)) * rstart rend resdn(i)
set res(i) rstart rend = %if(d1,resup(i),resdn(i))
end do i
*
* Compute the upper and lower branch covariance matrices, their Choleski
* factors and the inverse factor. (The inverse is for standardizing the
* residuals and the factor for mapping standardized residuals back to
* the true levels.)
*
vcv(matrix=vup)
# resup
compute sup =%decomp(vup)
compute siup=inv(sup)
vcv(matrix=vdn)
# resdn
compute sdn =%decomp(vdn)
compute sidn=inv(sdn)
*
* Compute the joint covariance matrix and its factor.
*
vcv(matrix=vsigma)
# res
compute s=%decomp(vsigma)
*
* Compute (jointly) standardized residuals. Since we only use these as a
* VECTOR across variables at T (never as individual series), we define
* it as a SERIES[VECTOR].
*
dec series[vect] stdu
gset stdu rstart rend = %if(d1,siup*%xt(resup,t),sidn*%xt(resdn,t))
*
* Save the original data since we'll overwrite it as part of the
* bootstrap.
*
dec vect[series] data(nvar)
do i=1,nvar
set data(i) = depvars(i){0}
end do i
*
* Define the switching formula. Because the bootstrap requires taking
* the standardized residuals and reflating them using regime-specific
* covariance matrices, the reflation part has to be incorporated into
* the formula. Thus TVARF(i) is (in effect) an identity given the (still
* to be created) bootres series.
*
do i=1,nvar
frml tvarf(i) depvars(i) = %if(thrfrml{d}>thresh,$
fitud(&i,1)+%dot(%xrow(sup,&i),bootres),$
fitud(&i,2)+%dot(%xrow(sdn,&i),bootres))
end do i
*
* Build the threshold var model using the switching equations and the
* definitional identity for the threshold variable.
*
*************************************************************************
*
* This is application specific. A different number of variables requires
* a different list of the tvarf's.
*
group tvar tvarf(1) tvarf(2) tvarf(3) tvarf(4) tvarf(5) thrfrml
*
*************************************************************************
*
* Figure out which set of entries we want. This is the simplest way
* to get this. d1{0} takes two values (0 or 1) but the order isn't known
* in advance since it will depend upon which value is seen first in the
* <<rstart>> to <<rend>> range. So we have to figure out which of values(1)
* and values(2) (and the corresponding entries(1) and entries(2)) is the
* one that we need, given the choice for <<upper>>.
*
panel(group=d1{0},id=values,identries=entries) d1 rstart rend
{
if upper==1.and.values(1)==1.or.upper==0.and.values(1)==0
compute rdates=entries(1),nrep=%size(rdates)
else
compute rdates=entries(2),nrep=%size(rdates)
}
*
* Set up target series. There is one for each combination of test sign
* and sizes on the shock, variable shocked and target variable.
*
declare real size sign
compute nexp=%size(shocksizes)*%size(shocksigns)
dec vect[rect[series]] girfs(nexp)
do i=1,nexp
dim girfs(i)(nvar,nvar)
clear(zeros) girfs(i)
end do i
*
* This is the working range for calculating the GIRF's
*
compute wstart=rstart,wend=rstart+horizon-1
*
* Outer loop is over the initial conditions, which walk through the data
* points in the desired regime. The residuals are bootstrapped by taking
* the standardized residuals (across the entire range), shuffling them,
* and reflating them based upon the current threshold value in the
* generated series.
*
infobox(action=define,lower=1,upper=nrep,progress) "Bootstrapping Across Initial Values"
do jrep=1,nrep
infobox(current=jrep)
*
* Copy observed data into depvar slots
*
compute basedate=rdates(jrep)
do i=1,nvar
set depvars(i) wstart-maxlag wstart-1 = data(i)(t-wstart+basedate)
end do i
*
* Loop over bootstrap replications
*
do krep=1,nkrep
*
* Generate the bootstrap shuffle
*
boot rentries wstart wend rstart rend
*
* Generate the base set of shocks by premultiplying the
* bootstrapped standardized shocks by the factor of the overall
* covariance matrix.
*
gset bootuv wstart wend = %if(%ranflip(.5),+1,-1)*stdu(rentries(t))
*
gset bootres wstart wend = bootuv
*
forecast(model=tvar,from=wstart,to=wend,results=base)
compute ifill=0
dofor sign = 1 -1
dofor size = 1 2
compute ifill=ifill+1
do jshock=1,nvar
*
* Patch over the component for which are computing the
* response with the selected size and sign.
*
compute bootres(wstart)=bootuv(wstart)
compute bootres(wstart)(jshock)=sign*size
*
forecast(model=tvar,from=wstart,to=wend,results=withshock)
do i=1,nvar
set girfs(ifill)(i,jshock) wstart wend = girfs(ifill)(i,jshock)+withshock(i)-base(i)
end do i
end do jshock
end do sign
end do size
end do krep
end do jrep
infobox(action=remove)
*
do k=1,nexp
do i=1,nvar
do j=1,nvar
set girfs(k)(i,j) wstart wend = girfs(k)(i,j)/(nkrep*nrep)
end do j
end do i
end do k
*
* Move the data back
*
do i=1,nvar
set depvars(i) = data(i)
end do i
*
*************************************************************************
*
dec vect[series] graphs(5)
compute klabels=||"+1 SD","+2 SD","-1 SD","-2 SD"||
*
spgraph(vfields=nvar,hfields=1,footer=$
"Figure 2. Response of Inflation Rate to Shocks, Conditional on Regime")
do j=1,nvar
do i=1,5
set graphs(i) 1 horizon = girfs(i)(1,j)(t+wstart-1)
end do i
graph(series=graphs,number=0,key=upright,klabels=klabels,$
header="Tight Regime: Response of Output",subheader="Shock to "+shortlabels(j))
end do j
spgraph(done)