' 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

