options crt double; name kleinfm2; ? ? Klein-I FIML - high precision benchmark ? See kleinfml.tsp for a much simpler version to reproduce the ? Calzolari-Panattoni benchmark coefficients and BHHH SEs to 4 digits ? with the regular FIML command (with 100 iterations). ? This benchmark is accurate to 11-12 digits, and also reproduces ? the Hessian and Rhat SEs. It uses the ML(HITER=N) command for ? efficient iteration. Coding the LogL for ML is a bit complicated, ? though. The Rhat SEs are not particularly useful, because they ? are usually not consistent for nonlinear FIML models, as explained ? by C-P. ? by Clint Cummins 12/1998 ? freq a; smpl 20,41; load; t = year-1931; ? ? Equations for Klein-I, with identities substituted in. ? Parameter names are the same as in Berndt's book. frml e1 cn - [a0 + a1*(w1+w2) + a2*(cn+i+g-tx-w1-w2) + a3*p(-1)]; frml e2 i - [b0 + b1*(cn+i+g-tx-w1-w2) + b2*p(-1) + b3*klag]; frml e3 w1 - [g0 + g1*(cn+i+g -w2) + g2*e(-1) + g3*t]; SMPL 21,41 ; PARAM A0-A3 B0-B3 G0-G3; ? ? Run 3SLS first to get starting values for FIML list ivs c p(-1) klag e(-1) t tx g w2; 3SLS(INST=ivs) e1-e3; ? ? FIML benchmark from Calzolari and Panattoni, Econometrica, ? May 1988, p.702. ? C-P TSP 4.4 11 digits ? A0 18.34 18.34327218344233 ? (SE BHHH) (12.88 ) (12.87924571526443 ) ? (SE Hess) ( 4.626) ( 4.625659300639248) ? (SE Rhat) ( 2.485) ( 2.485029947202310) ? A1 .8018 .8018443391624640 ? (.0842) (.0841774677755224) ? (.0445) (.0444945215159401) ? (.0359) (.0358931088318759) ? A2 -.2324 -.2323887662328997 ? (1.931 ) (1.930531045179115 ) ? ( .5806) ( .5806206557263284) ? ( .3120) ( .3119557287235352) ? A3 .3857 .3856730901020590 ? (1.083 ) (1.083100147822362 ) ? ( .3017) ( .3017459570226678) ? ( .2174) ( .2173571197498432) ? B0 27.26 27.26386576310186 ? (21.47 ) (21.47193027065058 ) ? ( 9.535) ( 9.534605265457136) ? ( 7.938) ( 7.937698265360596) ? B1 -.8010 -.8010060259538374 ? (2.334 ) (2.334490511501871 ) ? ( .8402) ( .8401881666046065) ? ( .4914) ( .4914212528868886) ? B2 1.052 1.051852141335067 ? (1.404 ) (1.403961123343230 ) ? ( .4244) ( .4243612037953698) ? ( .3525) ( .3524593026815088) ? B3 -.1481 -.1480990630401963 ? (.0992) (.0991710114662201) ? (.0468) (.0467958672191025) ? (.0299) (.0298546803887931) ? G0 5.794 5.794287580524262 ? (4.645) (4.645311665539802) ? (3.241) (3.240622249156222) ? (1.804) (1.804426514321850) ? G1 .2341 .2341176397831136 ? (.0953) (.0952877979733572) ? (.0950) (.0950131682279719) ? (.0488) (.0488180332688998) ? G2 .2847 .2846766802287325 ? (.0617) (.0617447495521507) ? (.0629) (.0628629674273714) ? (.0452) (.0452086728609125) ? G3 .2348 .2348346571073103 ? (.0776) (.0776415633099408) ? (.0565) (.0565279686364303) ? (.0345) (.0345002677103643) ? LogL* -83.3238 -83.32380966999975 ? (*LogL computed by ? TSP at the above ? coefficient values) ? ? We don't need to use the FIML command; it is faster to use ? ML at this point. ?FIML(endog=(cn,i,w1),maxit=100,tol=.0001) e1-e3; ? ? The ML code below is copied from fiml4.tsp ? ? instead of -.5*log(det(sigma)), make simplification: ? det(sigma) = 1/det(sigma-inverse) = 1/DL1*DL1*DL2*DL2*... ? log(det(sigma)) = -[2*log(DL1) + 2*log(DL2) + ...] ? = -2*ldetd (and -.5*(-2) = 1) ? or could use log(DL1*DL1) if there are negative problems ? set pi = 4*atan(1); set con = -.5*log(2*pi); frml fiml logl = ldetd + ldetjac -.5*epsie + con*neq; ? ? parameterize sigma inverse in terms of its LDL' factorization, ? so that its determinant is simplified ? s.. = sigma inverse element (lower triangle) ? d. = D element (diagonal matrix) ? DL. = square root of D element ? L.. = L element (unit triangular: ones on diagonal, zeros above diagonal) ? a.. = A element (A = LD) ? ? 2 x 2 case example: ? ? sigma inverse = s11 s21 = [ 1 0 ] * [ DL1*DL1 0 ] * [ 1 L21 ] ? s21 s22 = [ L21 1 ] * [ 0 DL2*DL2 ] * [ 0 1 ] ? ? These elements of sigma inverse are the same for higher order models; ? just add new rows to sigma inverse. ? frml d1 DL1*DL1; frml s11 d1; list es1 s11 d1; ? list of substitution equations list fp1 DL1 1; ? list of covariance parameters ? frml s21 s11*L21; frml d2 DL2*DL2; frml s22 s21*L21 + d2; list es2 s22 d2 s21 es1; list fp2 fp1 L21 DL2 1; ? frml s31 s11*L31; frml a32 d2*L32; frml s32 s31*L21 + a32; frml d3 DL3*DL3; frml s33 s31*L31 + a32*L32 + d3; list es3 s33 d3 s32 a32 s31 es2; list fp3 fp2 L31 L32 DL3 1; ? We can use abs() here to prevent numerical warnings of ? log of negative variance parameters. This will clean up ? the output a bit, but it won't speed it up any. In a hard ? coded implementation, these constraints would be checked directly. ?frml det3 ldetd = log(abs(DL1)) + log(abs(DL2)) + log(abs(DL3)); frml det3 ldetd = log(DL1) + log(DL2) + log(DL3); frml ese3 epsie = e1*(e1*s11 + 2*e2*s21 + 2*e3*s31) + e2*(e2*s22 + 2*e3*s32) + e3*e3*s33; ? ? log( abs (det Jacobian) ) for Klein I ? See Berndt exercise 10e8 to see how this expression for the determinant ? of the Jacobian was obtained. Or see the Jacobian equation in the Rhat ? section below. http://www.stanford.edu/~clint/berndt/10e8.tsp . ? The abs() is skipped because it is not needed this close to the optimum. ? frml Lj3 ldetjac = log( g1*(a2-a1+b1) + 1-a2-b1 ); set neq=3; param fp3; eqsub fiml det3 ese3 es3 Lj3 e3 e2 e1; ? ? Initialize covariance parameters from residual covariance ? matrix stored by 3SLS. ? inv @covu si; yldfac si dsig tsig; set DL1 = sqrt(dsig(1,1)); set DL2 = sqrt(dsig(2,2)); set DL3 = sqrt(dsig(3,3)); set L21 = tsig(1,2); set L31 = tsig(1,3); set L32 = tsig(2,3); ? ? The ML Hessian and BHHH SEs reproduce C-P options nwidth=23,signif=15; ml(hiter=n,hcov=nb,tol=1e-11) fiml; ? ? Run the FIML command to check its BHHH SEs ? (they reproduce C-P and ML). HCOV=G does not reproduce ? Rhat (the HCOV=G SEs are not normally consistent, so they ? are not normally printed). ? This will also store the residuals and residual covariance ? matrix at the optimum, for use in the Rhat calculations. ? FIML(endog=(cn,i,w1),hcov=bg) e1-e3; ? ? Form the Rhat matrix, following the details in the C-P ? discussion papers. Note: they follow Amemiya's notation, which ? assumes no cross-equation restrictions (correct for Klein-I). ? ?frml e1 cn - [a0 + a1*(w1+w2) + a2*(cn+i+g-tx-w1-w2) + a3*p(-1)]; ?frml e2 i - [b0 + b1*(cn+i+g-tx-w1-w2) + b2*p(-1) + b3*klag]; ?frml e3 w1 - [g0 + g1*(cn+i+g -w2) + g2*e(-1) + g3*t]; ? (-) derivatives with respect to parameters da1 = w1 + w2; da2 = cn + i + g - tx - w1 - w2; db1 = da2; dg1 = cn + i + g - w2; ? - G(i) in C-P notation mmake gb1 c da1 da2 p(-1); ? -G(1) mmake gb2 c db1 p(-1) klag; mmake gb3 c dg1 e(-1) t ; ? Jacobian mform(nrow=3,ncol=3,type=gen) j; set j(1,1) = 1-a2; set j(1,2) = -a2; set j(1,3) = a2-a1; set j(2,1) = -b1; set j(2,2) = 1-b1; set j(2,3) = b1; set j(3,1) = -g1; set j(3,2) = -g1; set j(3,3) = 1 ; mat ji = j"; ? inverse Jacobian read(nrow=4,ncol=3) gy1; ? dg1t/dyt' 0 0 0 0 0 -1 -1 -1 1 0 0 0 ; read(nrow=4,ncol=3) gy2; ? dg2t/dyt' 0 0 0 -1 -1 1 0 0 0 0 0 0 ; read(nrow=4,ncol=3) gy3; ? dg3t/dyt' 0 0 0 -1 -1 0 0 0 0 0 0 0 ; ? compute -Ghat(i) = F*(gy(i)*J")' - G(i) dot(value=i) 1-3; mat ghat. = @res*(gy.*ji)' + gb.; ? not - gb. because gb. = -G(i) enddot; mform(block) ghat = ghat1 ghat2 ghat3; mat sigi = @covu"; mat it = ident(@nob); mat sii = sigi # it; mat ghsi = ghat'sii; mat grad = ghsi*vec(@res); print grad; ? print gradient as check of calculation (should be zero) ?mat r = ghat'sii*ghat; mat r = sym(ghsi*ghat); mat ri = r"; ? reproduces C-P results for SEs from Rhat title 'Rhat SEs for FIML'; tstats(names=@rnms) @coef ri; end; noprint; freq a; smpl 20,41; load YEAR CN P W1 I KLAG E W2 G TX; 1920 39.8 12.7 28.8 2.7 180.1 44.9 2.2 4.6 3.4 1921 41.9 12.4 25.5 -0.2 182.8 45.6 2.7 6.6 7.7 1922 45.0 16.9 29.3 1.9 182.6 50.1 2.9 6.1 3.9 1923 49.2 18.4 34.1 5.2 184.5 57.2 2.9 5.7 4.7 1924 50.6 19.4 33.9 3.0 189.7 57.1 3.1 6.6 3.8 1925 52.6 20.1 35.4 5.1 192.7 61.0 3.2 6.5 5.5 1926 55.1 19.6 37.4 5.6 197.8 64.0 3.3 6.6 7.0 1927 56.2 19.8 37.9 4.2 203.4 64.4 3.6 7.6 6.7 1928 57.3 21.1 39.2 3.0 207.6 64.5 3.7 7.9 4.2 1929 57.8 21.7 41.3 5.1 210.6 67.0 4.0 8.1 4.0 1930 55.0 15.6 37.9 1.0 215.7 61.2 4.2 9.4 7.7 1931 50.9 11.4 34.5 -3.4 216.7 53.4 4.8 10.7 7.5 1932 45.6 7.0 29.0 -6.2 213.3 44.3 5.3 10.2 8.3 1933 46.5 11.2 28.5 -5.1 207.1 45.1 5.6 9.3 5.4 1934 48.7 12.3 30.6 -3.0 202.0 49.7 6.0 10.0 6.8 1935 51.3 14.0 33.2 -1.3 199.0 54.4 6.1 10.5 7.2 1936 57.7 17.6 36.8 2.1 197.7 62.7 7.4 10.3 8.3 1937 58.7 17.3 41.0 2.0 199.8 65.0 6.7 11.0 6.7 1938 57.5 15.3 38.2 -1.9 201.8 60.9 7.7 13.0 7.4 1939 61.6 19.0 41.6 1.3 199.9 69.5 7.8 14.4 8.9 1940 65.0 21.1 45.0 3.3 201.2 75.7 8.0 15.4 9.6 1941 69.7 23.5 53.3 4.9 204.5 88.4 8.5 22.3 11.6