* SHAZAM program for Full Information Maximum Likelihood (FIML).
* NOTE: This illustrates a method of programming an algorithm in SHAZAM.
* This is not necessarily recommended as a general way of doing FIML.
*
* This example was proposed by Professor Ray Byron of Bond University 
* in Australia. Although this method may not be the fastest 
* (or most appealing) way to compute FIML estimates it does illustrate 
* SHAZAM programming. 
* The example shows FIML estimation of Klein's model I. 
* The starting values for the iterative estimation are obtained from 
* 3SLS estimation (see the example in the chapter TWO-STAGE LEAST SQUARES 
* AND SYSTEMS OF EQUATIONS). The parameter estimates from the FIML 
* estimation are saved in the variable BFIML.
* -----------------------------------------------------------------------
* WG     GOVERNMENT WAGE BILL
* T      INDIRECT TAXES
* GE     GOVERNMENT EXPENDITURES
* TIME1  TIME TREND 
* PLAG   PROFITS LAGGED ONE PERIOD
* KLAG   STOCK OF CAPITAL AT END OF PREVIOUS PERIOD
* XLAG   PRIVATE PRODUCT LAGGED ONE PERIOD
* C      CONSUMPTION
* I      INVESTMENT
* WP     PRIVATE WAGE BILL
* PR     PROFITS
* WGWP   WAGES
* PP     PRIVATE PRODUCT
* -----------------------------------------------------------------------
SAMPLE 1 21
READ (klein.txt) WG T GE TIME1 PLAG KLAG XLAG C I WP PR WGWP PP 
* ===== KLEIN'S MODEL I - 3SLS =====
SYSTEM 3 WG T GE TIME1 PLAG KLAG XLAG / DN COEF=B3SLS
OLS C PLAG PR WGWP
OLS I PLAG KLAG PR
OLS WP XLAG TIME1 PP
* ===== KLEIN'S MODEL I - FIML =====
* Programmed by Ray Byron, Bond University, Queensland, Australia
* Note - this can be improved upon with
* two sided derivatives and better programming
*
* The model is  YB + XG = E
DIM B 7 7 G 7 7
* Exogenous variables
COPY PLAG KLAG XLAG TIME1 WG T GE  X
* Endogenous variables
GENR YY=C+I+GE-T
COPY C I WP YY PR WGWP PP  Y
* Number of observations
GEN1 NOB=21
* Number of stochastic equations
GEN1 NEQ=3
* Number of parameters
GEN1 NCOEF=9
GEN1 NM=NOB*NEQ
SET NODOECHO
* Specify the identities
*  (4) YY=C+I+GE-T
*  (5) PR=YY-WGWP
*  (6) WGWP=WP+WG
*  (7) PP=YY-WG+T
MATRIX B=-IDEN(7)
MATRIX G(6,4)=-1
MATRIX G(7,4)=1
MATRIX G(5,6)=1
MATRIX G(5,7)=-1
MATRIX G(6,7)=1
MATRIX B(1,4)=1
MATRIX B(2,4)=1
MATRIX B(4,5)=1
MATRIX B(6,5)=-1
MATRIX B(3,6)=1
MATRIX B(4,7)=1
* Express data in deviations about the mean
GENR ONE=1
MATRIX XM=ONE'X/NOB
MATRIX YM=ONE'Y/NOB
MATRIX X=X-ONE*XM
MATRIX Y=Y-ONE*YM
* Starting values
MATRIX BB=B3SLS
DIM F 3 1  DV NM NCOEF  E0 NM 1  E NM 1  V1 NOB NEQ  
* 
PROC ML_CALC
* Specify the stochastic equations
* (1) C = f(PLAG, PR, WGWP)     Consumption
* (2) I = f(PLAG, KLAG, PR)     Investment
* (3) WP= f(XLAG, TIME1, PP)    Private Wages
   MATRIX B(5,1)=BB(2)
   MATRIX B(6,1)=BB(3)
   MATRIX B(5,2)=BB(6)
   MATRIX B(7,3)=BB(9)
   MATRIX G(1,1)=BB(1)
   MATRIX G(1,2)=BB(4)
   MATRIX G(2,2)=BB(5)
   MATRIX G(3,3)=BB(7)
   MATRIX G(4,3)=BB(8)
   MATRIX P=-G*INV(B)
   MATRIX V=Y-X*P
   COPY V V1 / FROW=1;NOB TROW=1;NOB FCOL=1;NEQ TCOL=1;NEQ
   MATRIX VV=V1'V1/NOB
* DT is the value of the likelihood function
   MATRIX DT=DET(VV)
PROCEND
* 
EXEC ML_CALC
*
MATRIX FX1=DT
MATRIX BFIML=BB
MATRIX F(2)=DT
* set up matrix for solving quadratics later on
READ AA / ROWS=3 COLS=3
.01 -.1 1
 0   0  1
.01 .1  1
* loop for optimum searching
DO %=1,20
* routine to evaluate derivatives and get search direction
* vectorise v1
MATRIX E=VEC(V1)
MATRIX E0=E
MATRIX SIG=INV(VV)
* get dv/db numerically
DO #=1,NCOEF
* evaluate residuals in response to increments in coefficients
   MATRIX C=BFIML(#)*.01
   MATRIX BB(#)=BFIML(#)+C
   EXEC ML_CALC
   MATRIX E=VEC(V1)
   MATRIX DIF=(E-E0)/C
   COPY DIF DV / FROW=1;NM TROW=1;NM FCOL=1;1 TCOL=#;#
   MATRIX BB(#)=BFIML(#)
ENDO
* get weighting matrix and first derivatives for Gauss method
MATRIX DD=DV'(SIG@IDEN(NOB))*DV 
MATRIX D1=DV'(SIG@IDEN(NOB))*E
* get direction
MATRIX DIR=INV(DD)*D1
* test both sides of previous position
MATRIX BB=BFIML-.1*DIR
EXEC ML_CALC
MATRIX F(1)=DT
MATRIX BB=BFIML+.1*DIR
EXEC ML_CALC
MATRIX F(3)=DT
* solve quadratic
MATRIX ABC=INV(AA)*F
MATRIX LEN=-ABC(2)/(2*ABC(1))
MATRIX BB=BFIML+LEN*DIR
EXEC ML_CALC
GEN1 TEST=ABS((F(2)-DT)/F(2))
GEN1 F(2)=DT
MATRIX BFIML=BB
GEN1 DT1=1.0/DT
PRINT LEN DT1 BFIML
* Stopping criterion
ENDIF(TEST.LT..001)
ENDO
* Compare 3SLS and FIML
SAMPLE 1 NCOEF
PRINT B3SLS BFIML
STOP