Balke REStat 2000 |
These are replication files for the threshold VAR model in Balke(2000). This is a four variable VAR on quarterly U.S. data with a threshold break that separates into "tight credit" and "normal credit" regimes. The analysis allows for both the coefficients and covariance matrices to be different in the two regimes. A detailed description of the programs is included in the Structural Breaks and Switching Models e-course.
This includes two programs
tvar_estimate.rpf
Tests for a threshold break and estimates the model.
tvar_irf.rpf
Takes the estimated threshold value from tvar_estimate.rpf as given and computes and graphs non-linear impulse response functions using bootstrapping.
The author tries several different variables for the "credit" measure, which is both one of the variables in the model and used to create the threshold variable. These programs look at just one of the three: the spread between commercial paper rates and the Treasury bill rate. A high value for the spread creates the "tight" regime, as commercial customers have to pay a higher premium for credit.
The basic model (used in both programs) is
system(model=varmodel)
variables d1y d1p money credit
lags 1 to maxlag
det constant
end(system)
where DLY is the (annualized) growth rate in real GDP, DLP is the annualized inflation rate, MONEY is the Federal funds rate and CREDIT is the credit variable.
The threshold variable is a once-lagged moving average of credit variable. It has to be lagged (that is, predetermined) since it's generated from one of the dependent variables—without the delay, the model would be circular (threshold depends upon current value, current value depends upon threshold). The moving average is to smooth the signal out. The author used different moving averages for different credit measures; it's 2 for the spread variable. This generates this in a flexible fashion, using MALENGTH as an option for the moving average length.
filter(type=lagging,span=malength) credit / credthr
tvar_estimate.rpf
The optimal threshold value is determined by running systems regressions with different break values. The only places at which the model can break are at the observed values for the threshold series (CREDTHR), and the locations are also restricted to require at least a certain percentage of observations (here 15%, which is governed by the PI control value) over and above the minimum required to estimate the model. The estimates are done using the instruction SWEEP with the options VAR=HETERO (allowing for the covariances to different between regimes) and GROUP=CREDTHR{D}<=THRESH, which separates the entries into regimes based upon whether they are less than or equal to or are bigger than the test threshold value.
The estimated threshold value is 0.47, which is not what is shown in the paper—the result shown here is the best break for the full breaking model (both VAR coefficients and covariance matrices) that's used in the impulse response calculations, while the one in the paper used a model with some constraints. See: a graph of the likelihood ratio statistic vs the threshold.
As is typical in threshold models, the likelihood ratio has a non-standard distribution. The maximum likelihood ratio test statistic is 271.4 with 78 degrees of freedom (68 extra VAR coefficients and 10 in the covariance matrix), which is far out in the tails of the standard chi-squared, but the chi-squared isn't the proper distribution. Balke employs a fixed regressor bootstrap (Hansen(1996)); however, that doesn't apply to models where the covariance matrix changes with the regime. Our example does the test using fixed covariances instead. The results (subject to some simulation error) are:
Statistic P-value
Sup 164.9322 0.318
Avg 122.3309 0.122
Exp 79.3415 0.310
where the three statistics are three functionals on the likelihood ratio tests across break values from Andrews-Ploberger(1994): the sup (maximum), arithmetic average and exponential average. The p-values are the bootstrapped calculations, which are the percentage of the (fixed regressor) bootstrapped samples where that functional exceeded the one computed using the actual data. One would conclude that, if only the VAR coefficients are allowed to switch with the regime, there is no significant break in the data. Given that the sup with the covariance matrix changing is 271.4, compared to just 164.9322 without it, it looks like there much stronger evidence of a change in the covariance matrix than of a change in the coefficients.
tvar_irf.rpf
Impulse responses are the difference between forecasts of a system with and without certain added shocks. For linear models, this is independent of the initial conditions that govern the forecasts and the responses are linear in the added shocks (multiply the shocks by a constant, and the corresponding responses are multiplied by the same constant). Neither of those is true in the case of a non-linear model like this. To handle initial conditions, Balke does a straight average across all the observations is a particular regime. This example does the "upper" (or tight credit) regime. A different run would be needed to do the lower (or normal credit) regime.
comp upper =1
comp regimeDesc="Tight Regime"
He does 1 and 2 (standardized) unit shocks with each sign (since the responses don't "reflect" around 0 as they would with a linear model):
compute [vector] shocksizes=||+1.0,+2.0,-1.0,-2.0||
compute [vec[string]] shocksizel=||"+1 SD","+2 SD","-1 SD","-2 SD"||
This copies information about the threshold model from the other program (in particular, the optimal threshold value):
set credit = cpbill1
comp malength = 2
comp d = 1
comp thresh = 0.47000
This defines a VECTOR[INT] called DEPVARS with the set of dependent variables (which are used in defining the simulation model), with short and long version labels which will be used in the output. These will need to be changed if you apply to a different model.
compute depvars =||d1y,d1p,money,credit||
compute shortlabels=||"Output","Inflation","Money","Credit"||
compute longlabels =||"Output Growth","Inflation","Fed Funds","CP-Bill Spread"||
To do the estimation and testing in tvar_estimate.rpf, the threshold series (the lagged moving average of credit) can be treated as data. For out-of-sample calculations (forecasts and impulse responses) that's no longer the case—because it's a function of the dependent variable, it has to be calculated out-of-sample as the dependent variable changes. This sets up a (flexible) identity that creates the threshold series as the desired moving average of the CREDIT variable.
compute macoeffs=%fill(malength,1,1.0/malength)
equation(identity,coeffs=macoeffs) threqn credthr
# credit{0 to malength-1}
frml(equation=threqn,identity) thrfrml
This estimates the model under the two regimes, generating FITUP and FITLOFRML's for each equation for the upper and lower regime versions, RESUP and RESLO for the residuals in the upper and lower regimes and VUP and VLO for the covariance matrices in the upper and lower regimes.
estimate(smpl=dup,resids=resup,title="Upper Regime") rstart rend
do i=1,nvar
frml(equation=%modeleqn(basevar,i)) fitup(i)
end do i
compute vup=%sigma
*
estimate(smpl=dlo,resids=reslo,title="Lower Regime") rstart rend
do i=1,nvar
frml(equation=%modeleqn(basevar,i)) fitlo(i)
end do i
compute vlo=%sigma
Non-linear impulse responses require some form of simulation, as the forecasts themselves (and therefore the responses to added shocks) have no closed-form solutions and can only be done as averages across a "cloud" of simulated paths. Balke chooses to do parametric bootstrapping, using all the residuals across the two regimes. Since the covariance matrix differs between regimes, this requires jointly standardizing residuals in each regime to an identity covariance matrix, then "reflating" the reshuffled residuals to match the covariance matrix in the simulated regime. This defines the final formulas used for simulating each of the dependent variables given the simulated regime (controlled by THRFRML{D}>THRESH, which determines whether the upper or lower branch gets used), and the sampled standardized residuals. The &I notation is used for defining the FRML's in a loop.
do i=1,nvar
frml tvarf(i) depvars(i) = %if(thrfrml{d}>thresh,$
fitup(&i)+%dot(%xrow(sup,&i),bootres),$
fitlo(&i)+%dot(%xrow(slo,&i),bootres))
end do i
and this puts together the full model including the identity that defines THRFRML and also a separate one which defines a dummy for the current regime.
group tvar tvarf thrfrml upperfrml
Because the responses are computed by averaging across all observations which are classified into the regime that we want (here the upper/tight credit one), we need to generate a list of those observations. This is the simplest way to do that, relying on the PANEL instruction (which is also used for grouped data). In the end, RDATES will have a VECT[INTEGER] of the observations that we want.
panel(group=dup{0},id=values,identries=entries) dup rstart rend
{
if upper==1.and.values(1)==1.or.upper==0.and.values(1)==0
compute rdates=entries(1)
else
compute rdates=entries(2)
}
The bootstrap is run over the initial conditions on the outside and over the bootstrap shuffles on the inside. For a given bootstrap replication, this gets the base simulation, using a (sign) symmetrized reshuffling of the standardized residuals:
boot rentries wstart wend rstart rend
*
* Generate the base set of shocks through the simulation
* period by symmetrizing draws from the standardized
* residuals.
*
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)
In this case, there will now be 16 different impact shocks done one at a time: there are 4 variables (the JSHOCK loop) and 4 size/sign combinations for each (+/-1 and +/-2) (the K loop). This patches over the initial period shocks with the bootstrapped draws, then replaces the specific component with the desired value. Because these are all standardized, the results will be for +/-1 and +/-2 shocks in the Cholesky factor orthogonalized shock (as that's what we used earlier for factoring the covariance matrix).
compute bootres(wstart)=bootuv(wstart)
compute bootres(wstart)(jshock)=shocksizes(k)
The average of the difference between the forecasts with and without the added shock gives (an estimate of) the impulse response. How accurate that will be will depend upon the number of simulations. It will also depend upon how similar the initial conditions are to each other (or how insensitive the responses are to the initial conditions)—the dependence of the responses on initial conditions is something that can't be eliminated no matter how many draws you do.
The graphs of the impulse responses of output to the shocks in each variable are shown below. A common question is why there are no error bands. Note that simply computing the impulse responses themselves requires a fairly elaborate bootstrap—error bands would require another layer beyond that, and bootstrapping the entire threshold VAR is quite cumbersome. (This merely bootstraps taking the estimate TVAR as given).
Also shown is an estimated probability of being in the tight regime given shocks to the variables.
Graphs
From tvar_estimate.rpf, this shows the likelihood ratio statistic vs the threshold value.
.png)
From tvar_irf.rpf, this shows the responses of output growth to standardized shocks of various sizes and signs. Because of the Cholesky ordering, these will all be zero at step=0 except for output on output.
.png)
From tvar_irf.rpf, these show the simulated probability of being in the tight regime after shocks to each of the variables. (This only shows +/- 2 size shocks because the graphs would be too cluttered with all four combinations). Since this is conditional on starting in the tight regime, all have probability 1 at step=0.
.png)
Copyright © 2025 Thomas A. Doan