' 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

