The RATS Software Forum

Below find example code for using Uhlig's (1997, Econometrica) method to estimate a BVAR with stochastic volatility. Comments welcome.

************
************
************ Program including estimator of Uhlig's (Econometrica 1997) BVAR with stochastic volatility
************
************ Program written by Taisuke Nakata from the Federal Reserve Bank of Kansas City, with
************ edits by Todd Clark, based on Matlab code kindly provided by Harald Uhlig (for Uhlig estimator),
************ as well as Gibbs sampling BVAR code written by Tom Doan.
************
************ Using Uhlig's model and several others, the program generates forecasts of GDP growth, core PCE
************ inflation, and the federal funds rate. The VARs are specified in terms of GDP growth, core
************ inflation less an exponentially smoothed trend, and the funds rate less trend inflation
************ (such detrending makes the model similar to those considered in Kozicki and Tinsley (2001 JME, 2001 JEDC)).
************ In addition to displaying forecasts for 2007 and 2008, the program plots densities of forecasts from
************ the Uhlig model and constant volatility BVARs. The plots are saved to .eps files.
************
************ Note that, in contrast to the specifications used in Uhlig (1997), Doan, Litterman, and Sims (1984), among others,
************ the VARs are specified in terms of GDP growth, detrended inflation, and a real interest rate, rather than
************ log levels of GDP, the price index, and the nominal interest rate.
************
************ The program uses data on real GDP, total and core PCE price indexes, and the funds rate stored in a file data_q.xls
************
************ 02/09/2007
************
comp fixlags = 4 ;* fixed lag order to use in VARs

comp styr = 1947
cal styr 1 4
comp stsmpl = styr:1
comp endsmpl = 2006:4 ;* "endsmpl" is the last observation period.

all stsmpl endsmpl+10 ;* +10 to pad sample for forecasting
smpl stsmpl endsmpl

seed 2011+endsmpl ;* If you get two-hamped forecast distribution, assign a different value for "seed".

********************************
******************************** setting dimensions (# models, etc.)
********************************
comp nvar = 3 ;* number of variables in the VAR model (GDP, infl, ff rate)
comp nvar_real = 1 ;* number of real variables in the VAR model (GDP), ordered first in the system

comp nmeth = 6 ;* number of forecast methods to be considered
dec vec[string] methlabels(nmeth) ;* vector of strings to store short descriptions of each method

**** horizons (#s of periods)
comp [vec[int]] horz = ||1,2,4,8||
comp nhorz = %rows(horz)
comp maxh = horz(nhorz)

** note: the forecasts for 4 and 8 periods ahead are accumulated to obtain four-quarter rates of growth and inflation
** (but the quarterly level of the interest rate)

********************************
******************************** setting up and reading in data
********************************

open data data_q.xls
data(format=xls,org=col) / GDP PCEPI PCEXFEPI FFR
close

**** annualized log growth rates, stored in vector of series for use in VAR
dec vec[ser] y(nvar)
dofor i = GDP PCEXFEPI PCEPI
set(scratch) i = 400.*log(i{0}/i{1})
end do i

set infl = pcexfepi; set totalinfl = pcepi

**** storing data in vector of series for use in VAR
dec vec[ser] y(nvar)
comp j = 0
dofor i = GDP INFL FFR
comp j = j+1
set y(j) = i{0}
labels y(j)
# 'Y_'+%l(i)
end do i

**** differencing inflation and int rate, for use in benchmark MA models
dec vec[ser] dyser(nvar) ylevser(nvar)
dofor i = nvar_real+1 to nvar
set(scratch) dyser(i) = y(i){0} - y(i){1}
set(scratch) ylevser(i) = y(i){0}
end do i

**** estimate trend with exp. smoothing of inflation (yields a trend similar to the inflation expectations
**** series of Kozicki and Tinsley)
comp smparam = .07
inquire(series=infl) stpt_p *
inquire(series=totalinfl) stpt_totp *
clear ptrend
stats(noprint) totalinfl stpt_p-40 stpt_p-1
comp ptrend(stpt_p-1) = %mean

set ptrend stpt_p endsmpl = ptrend{1} + smparam*(infl - ptrend{1})
set rawptrend stpt_p endsmpl = ptrend ;* making a copy for use in the forecasting loop

set(scratch) y(nvar_real+1) = infl - ptrend{1}
set(scratch) y(nvar_real+2) = ffr - ptrend{1}

*print(dates) 1959:1 endsmpl infl ptrend

**** determining common sample to be used in estimation (after adjustment for lag order)
comp stpt = stsmpl
comp endpt = endsmpl
do i = 1,nvar
inquire(series=y(i)) stptser endptser
comp stpt = %imax(stptser,stpt)
comp endpt = %imin(endpt,endptser)
end do i

dis %datelabel(stpt) %datelabel(endpt)

********************************
******************************** stuff for Bayesian estimation:
********************************

******** prior means of first lag of dep variable in each equation (basic intention is to push VAR towards AR(1) models)
dec vec bvarprior(nvar)
comp bvarprior = %const(0.)
comp bvarprior(1) = .3 ;* GDP growth equation
comp bvarprior(nvar_real+1) = .7 ;* infl equation
comp bvarprior(nvar_real+2) = .9 ;* int rate equation

**** means for use in calibrating intercepts in Uhlig estimator
dec vec bvarmean(nvar)
comp bvarmean = %const(0.)
comp bvarmean(1) = 3.2 ;* GDP growth equation
comp bvarmean(nvar_real+1) = 0. ;* detrended infl equation
comp bvarmean(nvar_real+2) = 2.5 ;* (real) int rate equation

******** for conventional BVAR estimation

comp TIGHT = .2 ;* overall tightness of prior
comp OTHER = 0.5 ;* relative weight on other lags
comp decay = 1. ;* controls how the prior variance tightens as the lag order grows
comp lambda3 = decay ;* hyperparameter determining lag decay rate in full Minnesota implementation (the one using Gibbs sampling)

******* stuff for RATS' Gibbs sampling approach (taken from code by Tom Doan)
dec rect bover
dec vect bhat bmean bdiff ux
comp ndraws = 2000
comp burnindraws = ndraws/10
comp draws_burn = ndraws - burnindraws

******* stuff for Minnesota prior forecasts
declare rect sxx svtr swish betaols betadraw ranc
declare symm sigmad
dec vec univar(nvar)

********************************
******************************** stuff for Uhlig estimation
********************************
comp ndraws_uhlig = 4000
comp ncoef = nvar*fixlags + 1 ;*always include intercept

dec vec[rec] yser xser eser phiser ;* Needed to generate the data
dec vec[rec] nser invnser bhatser sser ;* Needed to calculate the posterior

comp niterate=21 ;* Following Uhlig's code.
dec vec[rec] bseq(niterate)
dec vec[vec] bseqv(niterate)
dec vec gradvec

smpl stsmpl endsmpl

dec vec udraw w_store(ndraws_uhlig)

***** stuff for Uhlig estimation
comp maxnobs = endsmpl-stsmpl+1
dim yser(maxnobs+maxh) xser(maxnobs+maxh) eser(maxnobs+1) phiser(maxnobs+1) nser(maxnobs+1) invnser(maxnobs+1) bhatser(maxnobs+1) sser(maxnobs+1)
do vtime=1,maxnobs
comp yser(vtime) = %zeros(nvar,1)
comp eser(vtime) = %zeros(nvar,1)
comp xser(vtime) = %zeros(ncoef,1)
comp nser(vtime) = %zeros(ncoef,ncoef)
comp invnser(vtime) = %zeros(ncoef,ncoef)
comp bhatser(vtime) = %zeros(nvar,ncoef)
comp sser(vtime) = %zeros(nvar,nvar)
end do vtime
comp v = 20. ;* recommended by Uhlig. (v=20 for quarterly data and v=60 for monthly data)
comp lambda = v/(v+1) ;* recommended by Uhlig.

** these zeta settings are meant to be in line with the conventional BVAR settings (and differ somewhat from Uhlig's suggestion)
comp zeta1 = 1./tight**2. ;* overall tightness of the prior
comp zeta2 = 1. ;* the increase in tightness for higher lags
comp zeta3 = 1./tight**2. ;* tightness of prior for constant
* (we impose an informative prior to prevent the intercept from moving around too much over time)

do vtime=1,maxnobs
comp xser(vtime)(1,1) = 1.
end do vtime

**** prior for coefs (same as in conventional BVAR)
do i=1,nvar
comp bhatser(1)(i,i+1) = bvarprior(i)
end do i
**** prior variance for intercept
comp nser(1)(1,1) = zeta3 ;* precision (the inverse of variance) of your prior for beta

**** prior mean for coefficients
dis ####.# bhatser(1)

********************************
******************************** stuff for forecasting
********************************

**** equation setup
dec vec[equations] vareqs(nvar) univareqs(nvar)
list mmm = 1 to nvar

dec vec[ser] forecasts(nvar) ;* vector of series used for temporary storage of quarterly forecasts
dec rec[real] fcy1(nmeth,nvar) fcy2(nmeth,nvar) fcy4(nmeth,nvar) fcy8(nmeth,nvar)
;* rectangular arrays for storing forecasts obtained with each method, one array for each horizon

** matrices for storing draws of Bayesian forecasts at the end of the sample, for density calculations
dec vec[rec] bfc1(nmeth) bfc2(nmeth) bfc4(nmeth) bfc8(nmeth) ;* separate vectors for each forecast horizon
do i = 1,nmeth
dim bfc1(i)(%imax(ndraws,ndraws_uhlig),nvar)
dim bfc2(i)(%imax(ndraws,ndraws_uhlig),nvar)
dim bfc4(i)(%imax(ndraws,ndraws_uhlig),nvar)
dim bfc8(i)(%imax(ndraws,ndraws_uhlig),nvar)
end do i

smpl stpt endsmpl+10
clear forecasts

comp stpt = stpt+fixlags ;* redefine starting point for sample estimate to reflect lags allowed

********************************
********************************
******************************** FORMING FORECASTS
********************************
********************************
comp rolln = 80 ;* size of sample used in rolling VAR estimates

comp time = endsmpl+1 ;* the forecast origin is endsmpl+1

dis '************************************************ forecast origin = ' %datelabel(time)
comp mod = 0
********************* defining common sample for VAR estimation
comp endpt = time-1
smpl stpt endpt
set ptrend stpt_p endsmpl = rawptrend{0}
*********************

********************* benchmark univariate AR(2) model for GDP growth, MA(1) models for changes in inflation, int rate
comp mod = mod+1
comp methlabels(mod) = 'univariate benchmark'
**
do i = 1,nvar_real
linreg(noprint,define=univareqs(i)) y(i)
# constant y(i){1 to 2}
ufore(equation=univareqs(i)) forecasts(i) endpt+1 endpt+maxh
end do i

comp rollstpt = %imax(stpt-fixlags+1,endpt-40+1) ;* start of estimation sample for rolling 10-year approach
dofor i = nvar_real+1 to nvar
boxjenk(ma=1,define=mamodel,noprint) dyser(i) rollstpt endpt
ufore(equation=mamodel) forecasts(i) endpt+1 endpt+maxh
accum forecasts(i) endpt+1 endpt+maxh accumfc
set forecasts(i) endpt+1 endpt+maxh = ylevser(i)(endpt) + accumfc{0}
end do i

*** forecasts for all variables
do i = 1,nvar
comp fcy1(mod,i) = forecasts(i)(time)
comp fcy2(mod,i) = forecasts(i)(time+1)
if i==nvar
{
comp fcy4(mod,i) = forecasts(i)(time+3)
comp fcy8(mod,i) = forecasts(i)(time+7)
}
else
{
comp fcy4(mod,i) = (forecasts(i)(time)+forecasts(i)(time+1)+forecasts(i)(time+2)+forecasts(i)(time+3))/4.
comp fcy8(mod,i) = (forecasts(i)(time+4)+forecasts(i)(time+5)+forecasts(i)(time+6)+forecasts(i)(time+7))/4.
}
end do i

********************* VAR
comp mod = mod+1
comp methlabels(mod) = 'VAR('+%string(fixlags)+')'
*

system(model=varmodel)
variables y
lags 1 to fixlags
kfset xxxols
det constant
end(system)
estimate(noprint)
forecast(model=varmodel,results=forecasts) * maxh endpt+1

*** storing stuff needed below in posterior simulations
compute sigma=%sigma
compute olssigma=%sigma
compute xx =inv(xxxols)
compute bhat=%vec(%modelgetcoeffs(varmodel))
compute sxx =%decomp(xxxols)
compute svtr =tr(%decomp(sigma))
compute betaols=%modelgetcoeffs(varmodel)
compute ncoef =%rows(sxx)
compute T = %nobs
compute wishdof= T - ncoef

*** obtaining forecasts of the level of infl and int rate from the detrended forecasts
set ptrend endpt+1 endpt+maxh = 0.
do ftime = endpt+1,endpt+maxh
comp forecasts(nvar-1)(ftime) = forecasts(nvar-1)(ftime) + ptrend(ftime-1)
comp forecasts(nvar)(ftime) = forecasts(nvar)(ftime) + ptrend(ftime-1)
comp ptrend(ftime) = ptrend(ftime-1) + smparam*(forecasts(nvar-1)(ftime) - ptrend(ftime-1))
end do ftime
*** forecasts for all variables
do i = 1,nvar
comp fcy1(mod,i) = forecasts(i)(time)
comp fcy2(mod,i) = forecasts(i)(time+1)
if i==nvar
{
comp fcy4(mod,i) = forecasts(i)(time+3)
comp fcy8(mod,i) = forecasts(i)(time+7)
}
else
{
comp fcy4(mod,i) = (forecasts(i)(time)+forecasts(i)(time+1)+forecasts(i)(time+2)+forecasts(i)(time+3))/4.
comp fcy8(mod,i) = (forecasts(i)(time+4)+forecasts(i)(time+5)+forecasts(i)(time+6)+forecasts(i)(time+7))/4.
}
end do i

********************* BVAR, Minnesota prior (using Gibbs sampling to get proper posterior mean of multi-step forecasts,
********************* as well as full posterior density)
********************* (note that I modified Doan's Gibbs sampling prior setup to follow the canned Minnesota settings)
comp mod = mod+1
comp methlabels(mod) = 'Minnesota BVAR('+%string(fixlags)+')'
do i = 1,nvar
comp fcy1(mod,i) = 0.
comp fcy2(mod,i) = 0.
comp fcy4(mod,i) = 0.
comp fcy8(mod,i) = 0.
end do i

*
* Create the prior mean and (square root of) precision. The olssee values are used to
* scale the precision. The prior mean is 1 for the first own lag in each equation, and
* 0 otherwise. The constant is given a mean of zero with a zero precision (that is, a flat
* prior).
*
comp lags = fixlags
compute ncoef=lags*nvar+1
comp minnmean = %zeros(ncoef,nvar)
comp minnprec = %zeros(ncoef,nvar)

** computing OLS univariate and VAR estimates, to get info needed in BVAR calcs
do i = 1,nvar
linreg(noprint) y(i)
# constant y(i){1 to lags}
compute univar(i)=sqrt(%rss/%nobs)
end do i

do i=1,nvar
do j=1,nvar
do l=1,lags
compute minnmean((j-1)*lags+l,i)=(i==j.and.l==1)*bvarprior(i)
compute minnprec((j-1)*lags+l,i)=univar(j)/univar(i)*%if(i==j,float(l)**lambda3/tight,float(l)**lambda3/(other*tight))
end do l
end do j
compute minnmean(ncoef,i)=0.0,minnprec(ncoef,i)=1./(1000.*univar(i))
end do i
*
* This works better with everything in vec form, rather than as an KxN matrix.
*
ewise minnprec(i,j)=minnprec(i,j)**2
compute hprior=%diag(%vec(minnprec))
compute bprior=%vec(minnmean)
*
* Except that there's one spot where the KxN form is more convenient, so we
* create bover and bdiff as overlayed copies of each other.
*
dim bmean(nvar*ncoef) bdiff(nvar*ncoef) ux(nvar*ncoef)
overlay bdiff(1) with bover(ncoef,nvar)
*
* The loop below does the Gibbs sampler. All I'm doing here is accumulating the means of
* the coefficients. However, the same basic code could easily be adapted for doing
* Monte Carlo integration for impulse response functions or forecasts.
*
do draws=1,ndraws
*
* Draw b given sigma. The draw for the coefficients is in bdraw as a vector.
*
compute hb=%kroneker(inv(sigma),xx)
compute hx=inv(hb+hprior)
compute bmean=hx*(hprior*bprior+hb*bhat)
compute fx=%decomp(hx)
compute bdraw=bmean+fx*(ux=%ran(1.0))
comp beta_draw = %vectorect(bdraw,ncoef)
*
* Draw sigma given b
*
compute bdiff=bdraw-bhat
compute wparm=olssigma+(1.0/%nobs)*%mqform(xx,bover)
compute sigma=%mqform(%nobs*inv(%ranwishart(nvar,%nobs-ncoef)),tr(%decomp(wparm)))
*
compute %modelsetcoeffs(varmodel,beta_draw)
simulate(model=varmodel,results=forecasts) * maxh endpt+1 sigma
*
if draws>burnindraws
{
*** obtaining forecasts of the level of infl and int rate from the detrended forecasts
set ptrend endpt+1 endpt+maxh = 0.
do ftime = endpt+1,endpt+maxh
comp forecasts(nvar-1)(ftime) = forecasts(nvar-1)(ftime) + ptrend(ftime-1)
comp forecasts(nvar)(ftime) = forecasts(nvar)(ftime) + ptrend(ftime-1)
comp ptrend(ftime) = ptrend(ftime-1) + smparam*(forecasts(nvar-1)(ftime) - ptrend(ftime-1))
end do ftime
*** forecasts for all variables
do i = 1,nvar
comp fcy1(mod,i) = fcy1(mod,i) + forecasts(i)(time)
comp fcy2(mod,i) = fcy2(mod,i) + forecasts(i)(time+1)
if i==nvar
{
comp fcy4(mod,i) = fcy4(mod,i) + forecasts(i)(time+3)
comp fcy8(mod,i) = fcy4(mod,i) + forecasts(i)(time+7)
}
else
{
comp fcy4(mod,i) = fcy4(mod,i) + (forecasts(i)(time)+forecasts(i)(time+1)+forecasts(i)(time+2)+forecasts(i)(time+3))/4.
comp fcy8(mod,i) = fcy4(mod,i) + (forecasts(i)(time+4)+forecasts(i)(time+5)+forecasts(i)(time+6)+forecasts(i)(time+7))/4.
}
end do i
*** storing for density calculations
do i = 1,nvar
comp bfc1(mod)(draws,i) = forecasts(i)(time)
comp bfc2(mod)(draws,i) = forecasts(i)(time+1)
if i==nvar
{
comp bfc4(mod)(draws,i) = forecasts(i)(time+3)
comp bfc8(mod)(draws,i) = forecasts(i)(time+7)
}
else
{
comp bfc4(mod)(draws,i) = (forecasts(i)(time)+forecasts(i)(time+1)+forecasts(i)(time+2)+forecasts(i)(time+3))/4.
comp bfc8(mod)(draws,i) = (forecasts(i)(time+4)+forecasts(i)(time+5)+forecasts(i)(time+6)+forecasts(i)(time+7))/4.
}
end do i

}
end do draws
*
do i = 1,nvar
comp fcy1(mod,i) = (1./draws_burn)*fcy1(mod,i)
comp fcy2(mod,i) = (1./draws_burn)*fcy2(mod,i)
comp fcy4(mod,i) = (1./draws_burn)*fcy4(mod,i)
comp fcy8(mod,i) = (1./draws_burn)*fcy8(mod,i)
end do i

********************* ROLLING VAR

*** first changing sample for rolling estimation
smpl %imax(stpt,endpt-rolln+1) endpt
***

comp mod = mod+1
comp methlabels(mod) = 'Rolling VAR('+%string(fixlags)+')'
*
****** ROLLING (OLS) estimates of VAR
system(model=varmodel)
variables y
lags 1 to fixlags
kfset xxxols
det constant
end(system)
estimate(print)
forecast(model=varmodel,results=forecasts) * maxh endpt+1

*** storing stuff needed below in posterior simulations
compute sigma=%sigma
compute olssigma=%sigma
compute xx =inv(xxxols)
compute bhat=%vec(%modelgetcoeffs(varmodel))
compute sxx =%decomp(xxxols)
compute svtr =tr(%decomp(sigma))
compute betaols=%modelgetcoeffs(varmodel)
compute ncoef =%rows(sxx)
compute T = %nobs
compute wishdof= T - ncoef

*** obtaining forecasts of the level of infl and int rate from the detrended forecasts
set ptrend endpt+1 endpt+maxh = 0.
do ftime = endpt+1,endpt+maxh
comp forecasts(nvar-1)(ftime) = forecasts(nvar-1)(ftime) + ptrend(ftime-1)
comp forecasts(nvar)(ftime) = forecasts(nvar)(ftime) + ptrend(ftime-1)
comp ptrend(ftime) = ptrend(ftime-1) + smparam*(forecasts(nvar-1)(ftime) - ptrend(ftime-1))
end do ftime
*** forecasts for all variables
do i = 1,nvar
comp fcy1(mod,i) = forecasts(i)(time)
comp fcy2(mod,i) = forecasts(i)(time+1)
if i==nvar
{
comp fcy4(mod,i) = forecasts(i)(time+3)
comp fcy8(mod,i) = forecasts(i)(time+7)
}
else
{
comp fcy4(mod,i) = (forecasts(i)(time)+forecasts(i)(time+1)+forecasts(i)(time+2)+forecasts(i)(time+3))/4.
comp fcy8(mod,i) = (forecasts(i)(time+4)+forecasts(i)(time+5)+forecasts(i)(time+6)+forecasts(i)(time+7))/4.
}
end do i

********************* ROLLING BVAR, Minnesota prior (using Gibbs sampling to get proper posterior mean of multi-step forecasts,
********************* as well as full posterior density)
comp mod = mod+1
comp methlabels(mod) = 'Rolling Minnesota BVAR('+%string(fixlags)+')'
do i = 1,nvar
comp fcy1(mod,i) = 0.
comp fcy2(mod,i) = 0.
comp fcy4(mod,i) = 0.
comp fcy8(mod,i) = 0.
end do i

*
* Create the prior mean and (square root of) precision. The olssee values are used to
* scale the precision. The prior mean is 1 for the first own lag in each equation, and
* 0 otherwise. The constant is given a mean of zero with a zero precision (that is, a flat
* prior).
*
comp lags = fixlags
compute ncoef=lags*nvar+1
comp minnmean = %zeros(ncoef,nvar)
comp minnprec = %zeros(ncoef,nvar)

** computing OLS univariate and VAR estimates, to get info needed in BVAR calcs
do i = 1,nvar
linreg(noprint) y(i)
# constant y(i){1 to lags}
compute univar(i)=sqrt(%rss/%nobs)
end do i

do i=1,nvar
do j=1,nvar
do l=1,lags
compute minnmean((j-1)*lags+l,i)=(i==j.and.l==1)*bvarprior(i)
compute minnprec((j-1)*lags+l,i)=univar(j)/univar(i)*%if(i==j,float(l)**lambda3/tight,float(l)**lambda3/(other*tight))
end do l
end do j
compute minnmean(ncoef,i)=0.0,minnprec(ncoef,i)=1./(1000.*univar(i))
end do i
*
* This works better with everything in vec form, rather than as an KxN matrix.
*
ewise minnprec(i,j)=minnprec(i,j)**2
compute hprior=%diag(%vec(minnprec))
compute bprior=%vec(minnmean)
*
* Except that there's one spot where the KxN form is more convenient, so we
* create bover and bdiff as overlayed copies of each other.
*
dim bmean(nvar*ncoef) bdiff(nvar*ncoef) ux(nvar*ncoef)
overlay bdiff(1) with bover(ncoef,nvar)
*
* The loop below does the Gibbs sampler. All I'm doing here is accumulating the means of
* the coefficients. However, the same basic code could easily be adapted for doing
* Monte Carlo integration for impulse response functions or forecasts.
*
do draws=1,ndraws
*
* Draw b given sigma. The draw for the coefficients is in bdraw as a vector.
*
compute hb=%kroneker(inv(sigma),xx)
compute hx=inv(hb+hprior)
compute bmean=hx*(hprior*bprior+hb*bhat)
compute fx=%decomp(hx)
compute bdraw=bmean+fx*(ux=%ran(1.0))
comp beta_draw = %vectorect(bdraw,ncoef)
*
* Draw sigma given b
*
compute bdiff=bdraw-bhat
compute wparm=olssigma+(1.0/%nobs)*%mqform(xx,bover)
compute sigma=%mqform(%nobs*inv(%ranwishart(nvar,%nobs-ncoef)),tr(%decomp(wparm)))
*
compute %modelsetcoeffs(varmodel,beta_draw)
simulate(model=varmodel,results=forecasts) * maxh endpt+1 sigma
*
if draws>burnindraws
{
*** obtaining forecasts of the level of infl and int rate from the detrended forecasts
set ptrend endpt+1 endpt+maxh = 0.
do ftime = endpt+1,endpt+maxh
comp forecasts(nvar-1)(ftime) = forecasts(nvar-1)(ftime) + ptrend(ftime-1)
comp forecasts(nvar)(ftime) = forecasts(nvar)(ftime) + ptrend(ftime-1)
comp ptrend(ftime) = ptrend(ftime-1) + smparam*(forecasts(nvar-1)(ftime) - ptrend(ftime-1))
end do ftime
*** forecasts for all variables
do i = 1,nvar
comp fcy1(mod,i) = fcy1(mod,i) + forecasts(i)(time)
comp fcy2(mod,i) = fcy2(mod,i) + forecasts(i)(time+1)
if i==nvar
{
comp fcy4(mod,i) = fcy4(mod,i) + forecasts(i)(time+3)
comp fcy8(mod,i) = fcy4(mod,i) + forecasts(i)(time+7)
}
else
{
comp fcy4(mod,i) = fcy4(mod,i) + (forecasts(i)(time)+forecasts(i)(time+1)+forecasts(i)(time+2)+forecasts(i)(time+3))/4.
comp fcy8(mod,i) = fcy4(mod,i) + (forecasts(i)(time+4)+forecasts(i)(time+5)+forecasts(i)(time+6)+forecasts(i)(time+7))/4.
}
end do i
*** storing for density calculations
do i = 1,nvar
comp bfc1(mod)(draws,i) = forecasts(i)(time)
comp bfc2(mod)(draws,i) = forecasts(i)(time+1)
if i==nvar
{
comp bfc4(mod)(draws,i) = forecasts(i)(time+3)
comp bfc8(mod)(draws,i) = forecasts(i)(time+7)
}
else
{
comp bfc4(mod)(draws,i) = (forecasts(i)(time)+forecasts(i)(time+1)+forecasts(i)(time+2)+forecasts(i)(time+3))/4.
comp bfc8(mod)(draws,i) = (forecasts(i)(time+4)+forecasts(i)(time+5)+forecasts(i)(time+6)+forecasts(i)(time+7))/4.
}
end do i
}
end do draws
*
do i = 1,nvar
comp fcy1(mod,i) = (1./draws_burn)*fcy1(mod,i)
comp fcy2(mod,i) = (1./draws_burn)*fcy2(mod,i)
comp fcy4(mod,i) = (1./draws_burn)*fcy4(mod,i)
comp fcy8(mod,i) = (1./draws_burn)*fcy8(mod,i)
end do i

********************** Uhlig BVAR
**********************
********************** General setup for Uhlig's procedure.
**********************
comp mod = mod+1
comp methlabels(mod) = 'Uhlig BVAR'
do i = 1,nvar
comp fcy1(mod,i) = 0.
comp fcy2(mod,i) = 0.
comp fcy4(mod,i) = 0.
comp fcy8(mod,i) = 0.
end do i

smpl stpt endpt
comp nobs = endpt - stpt + 1

do vtime=1,nobs
do i=1,nvar
comp yser(vtime)(i,1) = y(i)(stpt + vtime - 1)
do j=1,fixlags
comp test = y(i)(stpt+ vtime -j - 1)
comp xser(vtime)(1 + nvar*(j-1) + i,1) = test
end do j
end do i
end do vtime

**********************
********************** Specify the parameters for prior distribution of h(1)
**********************
** setting prior mean for intercept of growth, int. rate equation
dofor i = 1 nvar
stats(noprint) y(i) stpt endpt
comp bhatser(1)(i,1) = %mean*(1.-bvarprior(i))
end do i

do i=1,nvar
linreg(noprint) y(i) stpt endpt
# constant y(i){1 to fixlags}
comp sser(1)(i,i) = %rss/%nobs
end do i
*
do i=1,nvar
do j=1,fixlags
* precision (the inverse of variance) of your prior for beta
comp nser(1)(1+nvar*(j-1)+i,1+nvar*(j-1)+i) = (y(i)((stpt-fixlags))**2.)*zeta1*(j**zeta2)
end do j
end do i
comp invnser(1) = inv(nser(1))
*
do vtime=2,nobs+1
comp eser(vtime) = yser(vtime-1) - bhatser(vtime-1)*xser(vtime-1)
comp nser(vtime) = lambda*nser(vtime-1) + %outerxx(xser(vtime-1))
comp check1 = %mqform(%identity(1),(tr(xser(vtime-1))*invnser(vtime-1)))
comp check2 = %mqform(invnser(vtime-1),xser(vtime-1)) + lambda
comp invnser(vtime) = (invnser(vtime-1) - check1/check2(1,1))/lambda
comp bhatser(vtime) = (lambda*(bhatser(vtime-1)*nser(vtime-1)) + yser(vtime-1)*tr(xser(vtime-1)))*invnser(vtime)
comp [symm] xnx = %mqform(invnser(vtime),xser(vtime-1))
comp sser(vtime) = lambda*sser(vtime-1) + (lambda/v)*%outerxx(eser(vtime))*(1. - xnx(1,1))
end do vtime

********************
******************** We want to find an appropriate importance function.
********************
******* (1) Find "bmax", using modified Newton-Raphson method
*** (1-1) Calculate jstar(nvar*ncoef -by- nvar*ncoef)
comp jstar = %zeros(nvar*ncoef,nvar*ncoef)
do vtime=2,nobs
comp jstar = jstar - 0.5*%kroneker(nser(vtime),inv(v/lambda*sser(vtime)))
end do vtime
comp jstar = jstar - 0.5*(v+ncoef+1.)*%kroneker(nser(nobs +1),inv((v/lambda)*sser(nobs+1)))
comp invjstar = inv(jstar)

*** (1-2) Calculate a sequence of "B(n)" starting with B(1) = bhatser(T)
*
comp bseq(1) = bhatser(nobs+1)
comp test1 = %vectorect(%vec(bseq(1)),nvar)
*
do ni=2,niterate
comp bseq(ni) = %zeros(nvar,ncoef)
comp badj = %zeros(nvar,ncoef)
comp badjv = %zeros(nvar,ncoef)
comp gradvec= %zeros(nvar*ncoef,1)
*
do vtime=2,nobs+1
comp ae = inv( %mqform(nser(vtime),tr(bseq(ni-1)-bhatser(vtime))) + (v/lambda)*sser(vtime))
comp aevec = %zeros(1,nvar*nvar)
comp aevec = %vec(ae)
do i=1,nvar
do j=1,ncoef
comp entryij = (j-1)*nvar + i
comp eij = %zeros(nvar,ncoef)
comp eij(i,j)=1.
comp cij = eij*nser(vtime)*tr(bseq(ni-1)-bhatser(vtime)) + (bseq(ni-1)-bhatser(vtime))*nser(vtime)*tr(eij)
comp cijvec = %zeros(nvar*nvar,1)
comp cijvec = %vec(cij)
comp lltest = tr(aevec)*cijvec
if vtime<nobs+1
comp gradvec(entryij) = gradvec(entryij) - 0.5*lltest(1,1)
else
comp gradvec(entryij) = gradvec(entryij) - 0.5*(v+ncoef+1.)*lltest(1,1)
end do j
end do i
end do vtime
*
comp badjv = invjstar*gradvec
comp badj = %vectorect(badjv,nvar)
comp bseq(ni) = bseq(ni-1) - badj
end do ni
****** (1-3) We take bseq(ni) as the posterior mode.
comp bmax = bseq(ni)

dis ''
dis '********** Uhlig coef estimate = '
dis #####.### bmax
*
system(model=varmodel)
variables y
lags 1 to 4
det constant
kfset invzz
end(system)
estimate(noprint)
dis ''
dis '********** OLS coef estimate ='
dis #####.### %modelgetcoeffs(varmodel) ;*In comparing this with betadrawave, do not forget that the elements are orgnized differently.
;* They look roughly similar. Here, the intercept is in the last row.
;* Here, the matrix are blocked into different variables, and within each block,
;* you will see different lag length. In "bmax", the intercept is in the first column,
;* and the matrix are blocked into different lags, and futher into different variables.
comp secderiv=%zeros(nvar*ncoef,nvar*ncoef)

*** (2) Evaluate the Hessian of the marginal posterior (another parameter in the importance function)
do vtime=2,nobs+1
if vtime<nobs>= entryij
{
comp ekl = %zeros(nvar,ncoef)
comp ekl(k,l) = 1
comp ckl = ekl*nser(vtime)*tr(bmax-bhatser(vtime)) + (bmax-bhatser(vtime))*nser(vtime)*tr(ekl)
comp eneijkl = eij*nser(vtime)*tr(ekl) + ekl*nser(vtime)*tr(eij)
comp cklvec = %vec(ckl)
comp eneijklvec = %vec(eneijkl)
comp term1 = powerfactor*tr(aevec)*eneijklvec
comp term2 = powerfactor*tr(acavec)*cklvec
comp secderiv(entryij,entrykl) = secderiv(entryij,entrykl) + term1(1,1) - term2(1,1)
comp secderiv(entrykl,entryij) = secderiv(entryij,entrykl)
}
end do l
end do k
end do j
end do i
end do vtime

***** (3) Calculate adjustment factor for the weight and do remaining setup for importance sampling
comp vstar = fix(0.98*(nobs+v+ncoef-nvar*ncoef)) ;* the degrees of freedom in the importance function.
*
comp mp_add = 0.
comp mpbeta_bmax = 0.
do vtime=2,nobs+1
comp mp_add = 0.5*log( %det( %mqform(nser(vtime),tr(bmax-bhatser(vtime))) + (v/lambda)*sser(vtime) ) )
comp mpbeta_bmax = mpbeta_bmax + mp_add
end do vtime
comp mpbeta_bmax = mpbeta_bmax + 0.5*(v + ncoef)*log( %det( %mqform(nser(nobs+1),tr(bmax-bhatser(nobs+1))) + (v/lambda)*sser(nobs+1)) )
comp adj = 0.5*(vstar + nvar*ncoef)*log(1) - mpbeta_bmax;* the adjustment factor for weights
* ;* (includes the constant part of the marginal posterior of beta)
comp wbeta = 0.
comp beta = 0.
comp h = 0.
comp wh = 0.
comp wsum = 0.
comp invs = inv(secderiv)
comp invs_sqrt = %decomp(invs)
dim udraw(%rows(secderiv))
comp betadrawsum = %zeros(nvar,ncoef)
comp hdrawsum = %zeros(nvar,nvar)
comp wbetadrawsum = %zeros(nvar,ncoef)
comp whdrawsum = %zeros(nvar,nvar)

******* (4) Importance Sampling
do draws=1,ndraws_uhlig
********
******** Draw "betadraw" from t-distribution with mean=bmax, variance=jstar, and degree of freedom=vstar
********
comp betadrawv_demeaned = (0.5*vstar/%rangamma(0.5*vstar))*invs_sqrt*(udraw=%ran(1.0))
comp betadraw_demeaned = %vectorect(betadrawv_demeaned,nvar)
comp betadraw = bmax + betadraw_demeaned
********
******** Draw "hdraw" from Wishart-distribution with variance=sigma and shape parameter v + ncoef
********
comp invsigma = inv(lambda*%mqform(nser(nobs+1),tr(betadraw - bhatser(nobs+1))) + v*sser(nobs+1) )
comp cholinvsigma = tr(%decomp(invsigma))
comp hdraw = %mqform(%ranwishart(nvar,ncoef+v),cholinvsigma )
********
******** Generate forecast based on beta and H drawn (treating H as fixed over forecast horizon)
********
do ii=1,maxh
*
comp yser(nobs+ii) = %zeros(nvar,1)
comp xser(nobs+ii) = %zeros(ncoef,1)
*
comp xser(nobs+ii)(1,1) = 1.
do jj=2,1+nvar
comp xser(nobs+ii)(jj,1) = yser(nobs+ii-1)(jj-1,1)
end do jj
*
do jj=1+nvar+1,ncoef
comp xser(nobs+ii)(jj,1) = xser(nobs+ii-1)(jj-nvar,1)
end do jj
*
comp shocks = %ranmvnormal(%decomp(inv(hdraw)))
comp yser(nobs+ii) = betadraw*xser(nobs+ii) + shocks ;*%ranmvnormal(%identity(nvar))
end do ii
do ii = 1,nvar
set forecasts(ii) endpt+1 endpt+maxh = yser(nobs+t-endpt)(ii,1)
end do ii
*** obtaining forecasts of the level of infl and int rate from the detrended forecasts
set ptrend endpt+1 endpt+maxh = 0.
do ftime = endpt+1,endpt+maxh
comp forecasts(nvar-1)(ftime) = forecasts(nvar-1)(ftime) + ptrend(ftime-1)
comp forecasts(nvar)(ftime) = forecasts(nvar)(ftime) + ptrend(ftime-1)
comp ptrend(ftime) = ptrend(ftime-1) + smparam*(forecasts(nvar-1)(ftime) - ptrend(ftime-1))
end do ftime

***** Calculate weights
comp w_add = 0.
comp w_betadraw = 0.
do vtime=2,nobs+1
comp w_add = - 0.5*log( %det( %mqform(nser(vtime),tr(betadraw-bhatser(vtime))) + (v/lambda)*sser(vtime)) )
comp w_betadraw = w_betadraw + w_add
end do vtime
comp ttlg = (v + ncoef)*w_add
comp w_betadraw = w_betadraw + ttlg
comp term3 = %mqform(secderiv,%vec(betadraw_demeaned))
comp precimplg = 0.5*(vstar + nvar*ncoef)*log(1. + term3(1,1)/(vstar + nvar*ncoef))
comp w_betadraw = w_betadraw + precimplg
comp w_betadraw = w_betadraw - adj
comp w_betadraw = exp(w_betadraw)
comp wsum = wsum + w_betadraw
***** Weighted sum
comp wbetadrawsum = wbetadrawsum + betadraw*w_betadraw
comp whdrawsum = whdrawsum + hdraw*w_betadraw

*** forecasts for all variables
do i = 1,nvar
comp fcy1(mod,i) = fcy1(mod,i) + w_betadraw*forecasts(i)(time)
comp fcy2(mod,i) = fcy2(mod,i) + w_betadraw*forecasts(i)(time+1)
if i==nvar
{
comp fcy4(mod,i) = fcy4(mod,i) + w_betadraw*forecasts(i)(time+3)
comp fcy8(mod,i) = fcy4(mod,i) + w_betadraw*forecasts(i)(time+7)
}
else
{
comp fcy4(mod,i) = fcy4(mod,i) + w_betadraw*(forecasts(i)(time)+forecasts(i)(time+1)+forecasts(i)(time+2)+forecasts(i)(time+3))/4.
comp fcy8(mod,i) = fcy4(mod,i) + w_betadraw*(forecasts(i)(time+4)+forecasts(i)(time+5)+forecasts(i)(time+6)+forecasts(i)(time+7))/4.
}
end do i
*** storing for density calculations
do i = 1,nvar
comp bfc1(mod)(draws,i) = forecasts(i)(time)
comp bfc2(mod)(draws,i) = forecasts(i)(time+1)
if i==nvar
{
comp bfc4(mod)(draws,i) = forecasts(i)(time+3)
comp bfc8(mod)(draws,i) = forecasts(i)(time+7)
}
else
{
comp bfc4(mod)(draws,i) = (forecasts(i)(time)+forecasts(i)(time+1)+forecasts(i)(time+2)+forecasts(i)(time+3))/4.
comp bfc8(mod)(draws,i) = (forecasts(i)(time+4)+forecasts(i)(time+5)+forecasts(i)(time+6)+forecasts(i)(time+7))/4.
}
end do i
comp w_store(draws) = w_betadraw
end do draws
*
comp wbetadrawave = wbetadrawsum/wsum
comp whdrawave = whdrawsum/wsum
*
dis 'posterior mean of H = '
dis #####.### whdrawave
dis 'inverse of posterior mean of H = '
comp invH = inv(whdrawave)
dis #####.### invH

dis #####.### whdrawave
dis 'posterior mean of coefs = '
dis #####.### wbetadrawave
*
do i = 1,nvar
comp fcy1(mod,i) = (1./wsum)*fcy1(mod,i)
comp fcy2(mod,i) = (1./wsum)*fcy2(mod,i)
comp fcy4(mod,i) = (1./wsum)*fcy4(mod,i)
comp fcy8(mod,i) = (1./wsum)*fcy8(mod,i)
end do i


********************************
******************************** displaying forecasts as of the end of the sample
********************************
** 1-quarter ahead
do i = 1,nvar
if i==1
{
dis ''
dis '*********************** horizon = 1Q'
}
dis ''
dis '*************' %l(y(i))
dofor m = 1 to nmeth
dis methlabels(m) @48 ###.## fcy1(m,i)
end do m
end do i

** 2-quarter ahead
do i = 1,nvar
if i==1
{
dis ''
dis '*********************** horizon = 2Q'
}
dis ''
dis '*************' %l(y(i))
dofor m = 1 to nmeth
dis methlabels(m) @48 ###.## fcy2(m,i)
end do m
end do i

** 4-quarter ahead
do i = 1,nvar
if i==1
{
dis ''
dis '*********************** horizon = 1Y'
}
dis ''
dis '*************' %l(y(i))
dofor m = 1 to nmeth
dis methlabels(m) @48 ###.## fcy4(m,i)
end do m
end do i

** 8-quarter ahead
do i = 1,nvar
if i==1
{
dis ''
dis '*********************** horizon = 2Y'
}
dis ''
dis '*************' %l(y(i))
dofor m = 1 to nmeth
dis methlabels(m) @48 ###.## fcy8(m,i)
end do m
end do i


********************************
******************************** plotting posterior densities
********************************

grparm(portrait) header 18 subheader 18 keylabeling 16
env noshowgraphs gsave="BVARdens*.eps" gformat=landscape

set weightser 1 ndraws_uhlig = w_store(t)

*print 1 ndraws_uhlig weightser

***** 1-quarter ahead
spgraph(vfields=3,hfields=1,header='Posterior densities of 1Q-ahead forecasts')
do i = 1,nvar
if i==1
{
dis ''
dis '*************************************** horizon = 1Q'
}
dis ''
dis '*********************' %l(y(i))

dis '**** Minnesota prior, recursive:'
smpl burnindraws+1 ndraws
set statser = bfc1(3)(t,i)
stats(print,fractiles) statser
comp bw = (2.7*%nobs**(-.2))*%min(sqrt(%variance),((%fract75-%fract25)/1.34))
comp bwint = (%fract99-%fract01)/99.
density(bandwidth=bw,type=EPANECHNIKOV) statser / xvals3 yvals3

dis '**** Minnesota prior, rolling:'
smpl burnindraws+1 ndraws
set statser = bfc1(5)(t,i)
stats(print,fractiles) statser
comp bw = (2.7*%nobs**(-.2))*%min(sqrt(%variance),((%fract75-%fract25)/1.34))
comp bwint = (%fract99-%fract01)/99.
density(bandwidth=bw,type=EPANECHNIKOV) statser / xvals2 yvals2

dis '**** Uhlig'
smpl 1 ndraws_uhlig
set statser = bfc1(6)(t,i)
stats(print,fractiles) statser
comp bw = (2.7*%nobs**(-.2))*%min(sqrt(%variance),((%fract75-%fract25)/1.34))
comp bwint = (%fract99-%fract01)/99.
density(bandwidth=bw,type=EPANECHNIKOV,weight=weightser) statser / xvals yvals
*density(bandwidth=bw,type=EPANECHNIKOV) statser / xvals yvals

comp [string] header = %l(y(i))

scatter(style=line,nopatterns,header=header,key=below,klab=||'Uhlig','Minn-rolling','Minn-recursive'||) 3
# xvals yvals 1 100
# xvals2 yvals2 1 100
# xvals3 yvals3 1 100

end do i
spgraph(done)

***** 2-quarter ahead
spgraph(vfields=3,hfields=1,header='Posterior densities of 2Q-ahead forecasts')
do i = 1,nvar
if i==1
{
dis ''
dis '*************************************** horizon = 2Q'
}
dis ''
dis '*********************' %l(y(i))

dis '**** Minnesota prior, recursive:'
smpl burnindraws+1 ndraws
set statser = bfc2(3)(t,i)
stats(print,fractiles) statser
comp bw = (2.7*%nobs**(-.2))*%min(sqrt(%variance),((%fract75-%fract25)/1.34))
comp bwint = (%fract99-%fract01)/99.
density(bandwidth=bw,type=EPANECHNIKOV) statser / xvals3 yvals3

dis '**** Minnesota prior, rolling:'
smpl burnindraws+1 ndraws
set statser = bfc2(5)(t,i)
stats(print,fractiles) statser
comp bw = (2.7*%nobs**(-.2))*%min(sqrt(%variance),((%fract75-%fract25)/1.34))
comp bwint = (%fract99-%fract01)/99.
density(bandwidth=bw,type=EPANECHNIKOV) statser / xvals2 yvals2

dis '**** Uhlig'
smpl 1 ndraws_uhlig
set statser = bfc2(6)(t,i)
stats(print,fractiles) statser
comp bw = (2.7*%nobs**(-.2))*%min(sqrt(%variance),((%fract75-%fract25)/1.34))
comp bwint = (%fract99-%fract01)/99.
density(bandwidth=bw,type=EPANECHNIKOV,weight=weightser) statser / xvals yvals
*density(bandwidth=bw,type=EPANECHNIKOV) statser / xvals yvals

comp [string] header = %l(y(i))

scatter(style=line,nopatterns,header=header,key=below,klab=||'Uhlig','Minn-rolling','Minn-recursive'||) 3
# xvals yvals 1 100
# xvals2 yvals2 1 100
# xvals3 yvals3 1 100

end do i
spgraph(done)


***** 1-year ahead
spgraph(vfields=3,hfields=1,header='Posterior densities of 1Y-ahead forecasts')
do i = 1,nvar
if i==1
{
dis ''
dis '*************************************** horizon = 1Y'
}
dis ''
dis '*********************' %l(y(i))

dis '**** Minnesota prior, recursive:'
smpl burnindraws+1 ndraws
set statser = bfc4(3)(t,i)
stats(print,fractiles) statser
comp bw = (2.7*%nobs**(-.2))*%min(sqrt(%variance),((%fract75-%fract25)/1.34))
comp bwint = (%fract99-%fract01)/99.
density(bandwidth=bw,type=EPANECHNIKOV) statser / xvals3 yvals3

dis '**** Minnesota prior, rolling:'
smpl burnindraws+1 ndraws
set statser = bfc4(5)(t,i)
stats(print,fractiles) statser
comp bw = (2.7*%nobs**(-.2))*%min(sqrt(%variance),((%fract75-%fract25)/1.34))
comp bwint = (%fract99-%fract01)/99.
density(bandwidth=bw,type=EPANECHNIKOV) statser / xvals2 yvals2

dis '**** Uhlig'
smpl 1 ndraws_uhlig
set statser = bfc4(6)(t,i)
stats(print,fractiles) statser
comp bw = (2.7*%nobs**(-.2))*%min(sqrt(%variance),((%fract75-%fract25)/1.34))
comp bwint = (%fract99-%fract01)/99.
density(bandwidth=bw,type=EPANECHNIKOV,weight=weightser) statser / xvals yvals
*density(bandwidth=bw,type=EPANECHNIKOV) statser / xvals yvals

comp [string] header = %l(y(i))

scatter(style=line,nopatterns,header=header,key=below,klab=||'Uhlig','Minn-rolling','Minn-recursive'||) 3
# xvals yvals 1 100
# xvals2 yvals2 1 100
# xvals3 yvals3 1 100

end do i
spgraph(done)


***** 2-year ahead
spgraph(vfields=3,hfields=1,header='Posterior densities of 2Y-ahead forecasts')
do i = 1,nvar
if i==1
{
dis ''
dis '*************************************** horizon = 2Y'
}
dis ''
dis '*********************' %l(y(i))

dis '**** Minnesota prior, recursive:'
smpl burnindraws+1 ndraws
set statser = bfc8(3)(t,i)
stats(print,fractiles) statser
comp bw = (2.7*%nobs**(-.2))*%min(sqrt(%variance),((%fract75-%fract25)/1.34))
comp bwint = (%fract99-%fract01)/99.
density(bandwidth=bw,type=EPANECHNIKOV) statser / xvals3 yvals3

dis '**** Minnesota prior, rolling:'
smpl burnindraws+1 ndraws
set statser = bfc8(5)(t,i)
stats(print,fractiles) statser
comp bw = (2.7*%nobs**(-.2))*%min(sqrt(%variance),((%fract75-%fract25)/1.34))
comp bwint = (%fract99-%fract01)/99.
density(bandwidth=bw,type=EPANECHNIKOV) statser / xvals2 yvals2

dis '**** Uhlig'
smpl 1 ndraws_uhlig
set statser = bfc8(6)(t,i)
stats(print,fractiles) statser
comp bw = (2.7*%nobs**(-.2))*%min(sqrt(%variance),((%fract75-%fract25)/1.34))
comp bwint = (%fract99-%fract01)/99.
density(bandwidth=bw,type=EPANECHNIKOV,weight=weightser) statser / xvals yvals
*density(bandwidth=bw,type=EPANECHNIKOV) statser / xvals yvals

comp [string] header = %l(y(i))

scatter(style=line,nopatterns,header=header,key=below,klab=||'Uhlig','Minn-rolling','Minn-recursive'||) 3
# xvals yvals 1 100
# xvals2 yvals2 1 100
# xvals3 yvals3 1 100

end do i
spgraph(done)

end

I wanted to test your program but there is an error message in the middle of the execution:

## SX11. Identifier K is Not Recognizable. Incorrect Option Field or Parameter Order?
>>>>comp ekl(k,<<<<



It would be great if you could post an error free version of the program.
Alex

Sorry for the error. It seems that the process of pasting the text from my program file into the forum caused any lines spaced very far to the right to get deleted from the version posted on the forum. Sometime in the next few days I'll figure out how to post things correctly and thereby fix what is posted. In the meantime, though, if you email me at todd.e.clark@kc.frb.org, I can directly send you the original version that actually works.

Todd:

If you like, I can store a copy of it on our website and post a link here. Just email me a copy as an attachment if you want to do that.

Thanks,
Tom Maycock
support@estima.com

Thanks to Todd for emailing the files. You can download them using the following links.

Data file:
www.estima.com/procs/clark_uhlig/data_q.xls

Program file:
www.estima.com/procs/clark_uhlig/UhligBVAR.prg

From an earlier post in a different thread:

I am trying to apply Uhlig(1997) in a different context. I am a little unsure as to how to compute impulse responses at different points in time (say t1, t2 in a sample size T). The paper (and RATS codes by Todd) draws the coefficients and precision matrices at the last time period conditioned on data before that. If I were to compute impulse response at t1, does it mean that I need to restrict my estimation sample only until t1-1? And for t2, t2-1? This then looks like a rolling regression analysis. In other words, in Uhlig's scheme, we are not jointly drawing the entire history of parameters, but only parameters for a single time given data until that previous period. Is this correct?