Structural VARs

## Cholesky factorization as structural factorization

The following program illustrates the two forms of imposing the identifying restrictions in structural decompositions. For illustration purposes and to check that the restrictions are correctly imposed, we impose restrictions that replicate the Cholesky factorization. (If all you want is the standard impulse response based on the Cholesky factorization, there is no need to use the structrual decomposition.) The program checks whether the factorization matrix from the three methods match by computing the L-infinity norm (maximum of absolute column sums) of the matrix difference. Theoretically, the difference should be zero but because of floating errors the numerical difference won't be exactly zero.

    ' replicate cholesky factorization using SVAR
    ' for version 4
    ' 11/5/99 h
    ' last checked 10/27/2000 h

    ' include subroutine
    include sub_rmaxdiff.prg

    'change path to program path
    %path = @runpath
    cd "{%path}"

    ' create workfile
    workfile cholsvar q 1948:1 1979:3

    ' fetch data from database
    fetch(d=data_svar) rgnp rinv m1

    ' take log of each series
    series lgnp = log(rgnp)
    series linv = log(rinv)
    series lm1 = log(m1)

    ' estimate unrestricted VAR
    var var1.ls 1 4 lgnp linv lm1 @ c

    '-------------------------------------------------------------------
    ' method 1: short-run restrictions in text form
    '-------------------------------------------------------------------

    var1.cleartext(svar)
    var1.append(svar) @e1 = c(1)*@u1
    var1.append(svar) @e2 = -c(2)*@e1 + c(3)*@u2
    var1.append(svar) @e3 = -c(4)*@e1 - c(5)*@e2 + c(6)*@u3

    freeze(tab1) var1.svar(rtype=text,conv=1e-5)
    'show tab1

    ' store estimated A and B matrices from method 1
    matrix mata1 = var1.@svaramat
    matrix matb1 = var1.@svarbmat

    ' compute factorization matrix
    matrix fact1 = @inverse(mata1)*matb1

    '-------------------------------------------------------------------
    ' method 2: short-run restrictions in pattern matrix
    '-------------------------------------------------------------------

    ' create pattern matrices
    matrix(3,3) pata
    ' fill matrix in row major order
    pata.fill(by=r) 1,0,0, na,1,0, na,na,1

    matrix(3,3) patb
    patb.fill(by=r) na,0,0, 0,na,0, 0,0,na

    freeze(tab2) var1.svar(rtype=patsr,conv=1e-5,namea=pata,nameb=patb)
    'show tab2

    ' store estimated A and B matrices from method 2
    matrix mata2 = var1.@svaramat
    matrix matb2 = var1.@svarbmat

    ' compute factorization matrix
    matrix fact2 = @inverse(mata2)*matb2

    '-------------------------------------------------------------------
    ' method 3: built-in cholesky
    '-------------------------------------------------------------------

    ' do impulse response with default Cholesky ordering
    var1.impulse(10)

    ' store cholesky factor
    matrix fact3 = var1.@impfact

    '-------------------------------------------------------------------
    ' check whether results match
    '-------------------------------------------------------------------

    ' compute difference in factorization matrices
    scalar diff13
    call sub_rmaxdiff(fact1, fact3, diff13)
    scalar diff23
    call sub_rmaxdiff(fact2, fact3, diff23)

    ' display L-infinity norm of matrix difference in table
    table(2,2) check
    setcolwidth(check,1,20)
    check(1,1) = "method1 - method3:"
    check(1,2) = diff13
    check(2,1) = "method2 - method3:"
    check(2,2) = diff23

    show check

## Blanchard-Quah (1989) long-run restriction

The following program replicates the long-run restriction model by Blanchard-Quah (1989). (The data are not the same and the results won't exactly match.) The DY (output growth) and U (unemployment) series are already demeaned following the procedure by Blanchard-Quah (1989). The program estimates the factorization matrix by three different methods which should all yield the same result.

    ' Blanchard-Quah long-run restriction
    ' verify estimates using three methods
    ' for version 4
    ' 11/5/99 h
    ' last checked 10/27/2000 h

    ' include subroutine
    include sub_rmaxdiff.prg

    ' create workfile
    workfile blanquah q 1948:1 1987:4

    ' fetch data from database (already demeaned)
    fetch(d=data_svar) dy u

    ' estimate VAR (no constant)
    var var1.ls 1 8 dy u @

    '-------------------------------------------------------------------
    ' method 1: text form long-run restrictions
    '-------------------------------------------------------------------

    var1.cleartext(svar)
    var1.append(svar) @lr1(@u1)=0

    freeze(tab1) var1.svar(rtype=text,conv=1e-4)
    show tab1

    ' store estimated A and B matrices from method 1
    matrix mata1 = var1.@svaramat	' A should be identity
    matrix matb1 = var1.@svarbmat

    '-------------------------------------------------------------------
    ' method 2: long-run restrictions in short-run form
    ' *not recommended; only for checking*
    '-------------------------------------------------------------------

    ' get unit (cumulated) long-run response
    var1.impulse(imp=u)
    matrix clr = var1.@lrrsp

    var1.cleartext(svar)
    var1.append(svar) @e1 = c(1)*@u1 + c(2)*@u2
    var1.append(svar) @e2 = -clr(1,1)*c(1)/clr(1,2)*@u1 + c(4)*@u2

    freeze(tab2) var1.svar(rtype=text,conv=1e-5)
    show tab2

    ' store estimated A and B matrices from method 2
    matrix mata2 = var1.@svaramat	' A should be identity
    matrix matb2 = var1.@svarbmat

    '-------------------------------------------------------------------
    ' method 3: moment matching (used in Blanchard-Quah (1989) paper)
    ' *only works for just-identified models*
    '-------------------------------------------------------------------

    ' get residual covariance
    sym rcov = var1.@residcov
    ' and do cholesky factorization
    matrix chol = @cholesky(rcov)

    ' factorize unit long-run response
    matrix p = clr * chol

    ' factorized B matrix
    matrix(2,2) q
    ' impose sign normalization by taking positive root
    q(1,1) = @sqrt( 1/(1+p(1,1)*p(1,1)/p(1,2)/p(1,2)) )
    q(2,1) = -p(1,1)*q(1,1)/p(1,2)
    q(1,2) = -q(2,1)
    q(2,2) = q(1,1)

    ' get implied B matrix
    matrix matb3 = chol * q

    '-------------------------------------------------------------------
    ' check whether results match
    '-------------------------------------------------------------------

    ' difference in A matrices
    matrix eye2 = @identity(2)	' truth
    scalar adiff1
    call sub_rmaxdiff(mata1, eye2, adiff1)
    scalar adiff2
    call sub_rmaxdiff(mata2, eye2, adiff2)

    ' difference in B matrices
    scalar bdiff13
    call sub_rmaxdiff(matb1, matb3, bdiff13)
    scalar bdiff23
    call sub_rmaxdiff(matb2, matb3, bdiff23)

    ' display L-infinity norm of matrix difference in table
    table(2,3) check
    setcolwidth(check,1,10)
    setcolwidth(check,2,15)
    setcolwidth(check,3,15)

    check(1,2) = "A - eye2"
    check(1,3) = "B - method3"
    check(2,1) = "method1:"
    check(3,1) = "method2:"

    check(2,2) = adiff1
    check(3,2) = adiff2
    check(2,3) = bdiff13
    check(3,3) = bdiff23

    show check

## Short-run restriction impulse responses (Sims 1986)

The following program replicates the short-run restriction exercises in Sims (1986). (The data are not the same and the results do not match exactly.) Note that while Sims assumes a diagonal covariance matrix for the structural innovations, EViews assumes an identity covariance matrix. For this reason, we need to estimate the standard deviation of the structural shocks as elements of the B matrix.

    ' replicate Sims (1986) with similar data
    ' for version 4
    ' 1/12/2000 h
    ' last checked 10/27/2000 h

    'change path to program path
    %path = @runpath
    cd "{%path}"

    ' create workfile
    workfile sims1 q 1948:1 1979:3

    ' fetch data from database
    fetch(d=data_svar) rgnp rinv pidgnp m1 tb3sec unemp

    ' transform series
    series lgnp = log(rgnp)
    series linv = log(rinv)
    series lprc = log(pidgnp)
    series lm1 = log(m1)

    '----------------------------------------------------------------------
    ' replicate Chart 1 (p.11) *no Bayesian priors imposed*
    '----------------------------------------------------------------------

    var var1.ls 1 4 lgnp linv lprc lm1 unemp tb3sec @ c
    freeze(chart1) var1.impulse(32,m,imp=chol,se=a)
    show chart1

    ' plot only some of money innovation responses
    freeze(chart1a) var1.impulse(32,m,imp=chol,se=a) lgnp lprc tb3sec @ lm1

    '----------------------------------------------------------------------
    ' replicate Chart 2 (p.13)
    '----------------------------------------------------------------------

    ' reorder VAR to confirm to Sims' (1986) ordering
    var var1.ls 1 4 tb3sec lm1 lgnp lprc unemp linv @ c

    ' 1st identification restrictions (p.12 (12)-(17))
    coef(6) b
    var1.append(svar) @e1 = c(1)*@e2 + b(1)*@u1
    var1.append(svar) @e2 = c(2)*@e3 + c(3)*@e4 + c(4)*@e1 + b(2)*@u2
    var1.append(svar) @e3 = c(5)*@e1 + c(6)*@e6 + b(4)*@u4
    var1.append(svar) @e4 = c(7)*@e1 + c(8)*@e3 + c(9)*@e6 + b(5)*@u5
    var1.append(svar) @e5 = c(10)*@e1 + c(11)*@e3 + c(12)*@e6 + c(13)*@e4 + b(6)*@u6
    var1.append(svar) @e6 = b(3)*@u3

    ' estimate structural factorization
    freeze(svar1) var1.svar(rtype=text,conv=1e-5)
    show svar1

    ' replicate Chart 2 (p.13)
    freeze(chart2) var1.impulse(32,m,imp=struct,se=a) lgnp linv lprc lm1 unemp tb3sec @ 1 2 4 5 6 3
    show chart2

    ' plot only some of money supply responses
    freeze(chart2a) var1.impulse(32,m,imp=struct,se=a) lgnp lprc tb3sec @ 1

    '----------------------------------------------------------------------
    ' replicate Chart 3 (p.13)
    '----------------------------------------------------------------------

    ' clear previous restrictions
    var1.cleartext(svar)

    ' 2nd identification restrictions (p.14 (18)-(23))
    var1.append(svar) @e1 = c(1)*@e2 + b(1)*@u1
    var1.append(svar) @e2 = c(2)*@e3 + c(3)*@e6 + c(4)*@e4 + c(5)*@e1 + b(2)*@u2
    var1.append(svar) @e3 = c(6)*@e1 + c(7)*@e6 + b(3)*@u3
    var1.append(svar) @e4 = c(8)*@e1 + c(9)*@e3 + c(10)*@e2 + b(5)*@u5
    var1.append(svar) @e5 = c(11)*@e1 + c(12)*@e3 + c(13)*@e4 + c(14)*@e6 + b(6)*@u6
    var1.append(svar) @e6 = b(4)*@u4

    ' estimate structural factorization
    ' default starting value does not converge
    rndseed 123456
    freeze(svar2) var1.svar(rtype=text,conv=1e-5,f0=u)
    show svar2

    ' replicate Chart 3 (p.13)
    freeze(chart3) var1.impulse(32,m,imp=struct,se=a) lgnp linv lprc lm1 unemp tb3sec @ 1 2 3 5 6 4
    show chart3

    ' plot only some of money supply responses
    freeze(chart3a) var1.impulse(32,m,imp=struct,se=a) lgnp lprc tb3sec @ 1

## Long-run restriction impulse responses (Blanchard-Quah 1989)

The following program replicates the long-run restriction model by Blanchard-Quah (1989). (The data are not the same and the results do not exactly match.) The DY (output growth) and U (unemployment) series are already demeaned following the procedure in Blanchard-Quah (1989). In order to plot the accumulated response and the level response in one graph, we first store the impulse responses in separate matrices, extract and place the relevant columns into a third matrix and use the line view of that matrix.

    ' replicate impulse response graph
    ' from long-run restriction (Blanchard-Quah 1989)
    ' for version 4
    ' 11/12/99 h
    ' last checked 10/27/2000 h

    ' include subroutine
    include sub_rmaxdiff.prg

    ' create workfile
    workfile blanquah q 1948:1 1987:4

    ' fetch data from database (already demeaned)
    fetch(d=data_svar) dy u

    ' change to percentage
    dy = 100*dy

    ' estimate VAR (no constant)
    var var1.ls 1 8 dy u @

    ' impose long-run restrictions
    ' @u1 = aggregate demand shock
    ' @u2 = aggregate supply shock
    var1.cleartext(svar)
    var1.append(svar) @lr1(@u1)=0

    ' estimate factorization
    freeze(tab1) var1.svar(rtype=text,conv=1e-5)

    ' set response horizon
    !hrz = 40

    ' replicate Figures 3-4 (p.662)
    ' accumulate to get level response of output
    freeze(fig3) var1.impulse(!hrz,imp=struct,se=a,a,matbys=rsp_acc) dy @ 1
    freeze(fig4) var1.impulse(!hrz,imp=struct,se=a,a) dy @ 2
    freeze(fig34) fig3 fig4
    show fig34

    ' replicate Figures 5-6
    freeze(fig5) var1.impulse(!hrz,imp=struct,se=a,matbys=rsp_lvl) u @ 1
    freeze(fig6) var1.impulse(!hrz,imp=struct,se=a) u @ 2
    freeze(fig56) fig5 fig6
    show fig56

    ' replicate Figure 1
    matrix(!hrz,2) mtmp
    ' get accumulated output response
    vector v = @columnextract(rsp_acc,1)
    colplace(mtmp,v,1)
    ' get level unemployment response
    v = @columnextract(rsp_lvl,2)
    colplace(mtmp,v,2)

    freeze(fig1) mtmp.line
    fig1.draw(line,left) 0
    fig1.elem(1) legend(Output response to demand)
    fig1.elem(2) legend(Unemployment response to demand)
    show fig1

    ' replicate Figure 2
    ' get accumulated output response
    vector v = @columnextract(rsp_acc,3)
    colplace(mtmp,v,1)
    ' get level unemployment response
    v = @columnextract(rsp_lvl,4)
    colplace(mtmp,v,2)

    freeze(fig2) mtmp.line
    fig2.draw(line,left) 0
    fig2.elem(1) legend(Output response to supply)
    fig2.elem(2) legend(Unemployment response to supply)
    show fig2

## Bootstrap impulse response confidence bounds for structural VARs with long-run restrictions

Currently EViews does not provide standard error bounds for impulse responses from structural impulses identified by long-run restrictions. The following program illustrates how to obtain bootstrap impulse response confidence bounds for the Blanchard-Quah (1989) model. The bootstrap is performed as follows:

1. Resample from the VAR residuals with replacement.
2. Obtain bootstrap data using the estimated VAR coefficients and the bootstrap residuals.
3. Estimate the VAR and structural factorization using the bootstrap data from step 2.
4. Store the impulse responses using estimates from step 3.
5. Repeat steps 1-4 many times.

Step 2 is implemented by specifying the bootstrap residuals as (level shift) add factors in model solve. Note also that estimation of the structural factorization at each bootstrap step may fail to converge; in such cases the code only stores results if convergence was achieved. The code below does 100 bootstrap replications; you should probably increase this number to get more reliable estimates but be forewarned that this make take a very long time to execute.

    ' bootstrap structural VAR impulse responses
    ' from long-run restriction (Blanchard-Quah 1989)
    ' for version 4
    ' 11/28/2000 h

    include sub_bootirf.prg

    ' create workfile
    workfile blanquah q 1948:1 1987:4

    ' fetch data from database (already demeaned)
    fetch(d="i:\hiroyuki\example files\var\data_svar") dy u

    ' change to percentage
    dy = 100*dy

    ' estimate VAR (no constant)
    var var1.ls 1 8 dy u @

    ' impose long-run restrictions
    ' @u1 = aggregate demand shock
    ' @u2 = aggregate supply shock
    var1.cleartext(svar)
    var1.append(svar) @lr1(@u1)=0

    ' estimate factorization
    var1.svar(rtype=text,conv=1e-5)

    ' set response horizon
    !hrz = 40

    ' store level and accumulated responses
    do var1.impulse(!hrz,imp=struct,matbys=rsp_lvl) dy u @ 1
    do var1.impulse(!hrz,imp=struct,a,matbys=rsp_acc) dy u @ 1

    ' get residuals to bootstrap
    var1.makeresid(n=gres) res_dy res_u

    ' make model to generate bootstrap data
    var1.makemodel(mod1)
    ' assign add factors for bootstrap residuals
    mod1.addassign(i) @all

    'set random number generator
    rndseed(type=mt) 1234567
    'set number of bootstrap replications
    !reps = 100
    'declare storage matrices
    for !i=1 to @columns(rsp_lvl)
    	matrix(!hrz,!reps) brsp_lvl{!i}
    	matrix(!hrz,!reps) brsp_acc{!i}
    next
    'temporary working vector
    vector vtmp
    'count number of cases that converged
    scalar cnt = 0

    'bootstrap loop
    for !i=1 to !reps
    	'bootstrap residuals
    	smpl @all if res_dy<>na
    	gres.resample dy_a u_a
    	'generate bootstrap data
    	mod1.solve
    	'estimate VAR (no constant) with bootstrap data
    	smpl @all
    	var var2.ls 1 8 dy_0 u_0 @
    	'impose long-run restrictions
    	var2.cleartext(svar)
    	var2.append(svar) @lr1(@u1)=0
    	'estimate factorization
    	var2.svar(rtype=text,conv=1e-5,maxiter=100,nostop)
    	'store results only if converged
    	if (var2.@svarcvgtype=0) then
    		'increment counter
    		cnt = cnt + 1
    		'get level and accumulated bootstrap responses
    		do var2.impulse(!hrz,imp=struct,matbys=tmp_lvl) dy_0 u_0 @ 1
    		do var2.impulse(!hrz,imp=struct,a,matbys=tmp_acc) dy_0 u_0 @ 1
    		'and store in appropriate position
    		for !j=1 to @columns(rsp_lvl)
    			vtmp = @columnextract(tmp_lvl,!j)
    			colplace(brsp_lvl{!j}, vtmp, cnt)
    			vtmp = @columnextract(tmp_acc,!j)
    			colplace(brsp_acc{!j}, vtmp, cnt)
    		next
    		delete tmp_lvl tmp_acc
    	endif
    next

    'reshape results matrix if there were non-converged cases
    if (cnt <> !reps) then
    	for !i=1 to @columns(rsp_lvl)
    		brsp_lvl{!i} = @subextract(brsp_lvl{!i}, 1, 1, !hrz, cnt)
    		brsp_acc{!i} = @subextract(brsp_acc{!i}, 1, 1, !hrz, cnt)
    	next
    	'display warning in table
    	table warning
    	%msg = @str(!reps-cnt) + " case(s) out of " + @str(!reps) + " replications failed to converge"
    	warning(1,1) = %msg
    	show warning
    endif

    'find bootstrap 66 percentiles for each horizon
    matrix(!hrz,@columns(rsp_lvl)) lvl_upp
    matrix(!hrz,@columns(rsp_lvl)) lvl_low
    matrix(!hrz,@columns(rsp_lvl)) acc_upp
    matrix(!hrz,@columns(rsp_lvl)) acc_low

    rowvector rtmp
    for !i=1 to !hrz
    	for !j=1 to @columns(rsp_lvl)
    		rtmp = @rowextract(brsp_lvl{!j},!i)
    		lvl_upp(!i,!j) = @quantile(rtmp,0.833)
    		lvl_low(!i,!j) = @quantile(rtmp,0.167)
    		rtmp = @rowextract(brsp_acc{!j},!i)
    		acc_upp(!i,!j) = @quantile(rtmp,0.833)
    		acc_low(!i,!j) = @quantile(rtmp,0.167)
    	next
    next

    'plot impulse responses with bootstrap error bounds
    %gname = "fig_lvl"
    call sub_bootirf(rsp_lvl, lvl_upp, lvl_low, %gname)
    ' add zero line to all graphs
    {%gname}.draw(line,left,rgb(180,180,180)) 0
    ' add title to graph
    {%gname}.addtext(0.1,-0.2) Bootstrap 66% Confidence Bounds (Levels)
    ' realign
    {%gname}.align(2,1,0.5)
    show {%gname}

    %gname = "fig_acc"
    call sub_bootirf(rsp_acc, acc_upp, acc_low, %gname)
    ' add zero line to all graphs
    {%gname}.draw(line,left,rgb(180,180,180)) 0
    ' add title to graph
    {%gname}.addtext(0.1,-0.2) Bootstrap 66% Confidence Bounds (Accumulated Responses)
    ' realign
    {%gname}.align(2,1,0.5)
    show {%gname}