R version 2.10.1 (2009-12-14) Copyright (C) 2009 The R Foundation for Statistical Computing ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > ## ---------------------------------------------------------------------------- > ## Zellner (1962) > ## > ## data > ## availability: none > ## firms: General Electric, Westinghouse > ## errors: (none) > ## > ## analysis > ## result: OLS exact at precision given > ## only minor deviations at last significant digit for some variances > ## SUR estimates not reproducible due to algebraic error on p. 359 > ## ---------------------------------------------------------------------------- > > source("start.R") Loading required package: kinship Loading required package: survival Loading required package: splines Loading required package: nlme Loading required package: lattice [1] "kinship is loaded" Loading required package: Formula Loading required package: MASS Loading required package: sandwich Loading required package: zoo Loading required package: Matrix Loading required package: car Loading required package: lmtest > > gr <- subset(Grunfeld, firm %in% c("General Electric", "Westinghouse")) > gr$firm <- factor(gr$firm) > pgr <- plm.data(gr, c("firm", "year")) > > ## Two-Stage Aitken Method (Table 1, p. 360) > ## ?? > > ## OLS estimation (Table 1, p. 360) > fmGE <- lm(invest ~ value + capital, data = gr, subset = firm == "General Electric") > fmWH <- lm(invest ~ value + capital, data = gr, subset = firm == "Westinghouse") > diag(vcov(fmGE)) (Intercept) value capital 9.843435e+02 2.423036e-04 6.606999e-04 > diag(vcov(fmWH)) (Intercept) value capital 6.424486e+01 2.466942e-04 3.147095e-03 > ## alternatively: > ## fm <- systemfit(invest ~ value + capital, method = "OLS", data = pgr) > > > ## computations leading to two-stage estimates > ## p. 358 > mmGE <- as.matrix(cbind(model.frame(fmGE), 1)) > mmWH <- as.matrix(cbind(model.frame(fmWH), 1)) > > ## moment matrices > crossprod(mmGE) invest value capital 1 invest 254113.5 4093308.3 1005863.5 2045.8 value 4093308.3 78628914.2 15769824.1 38826.5 capital 1005863.5 15769824.1 4395946.8 8003.2 1 2045.8 38826.5 8003.2 20.0 > crossprod(mmGE, mmWH) invest value capital 1 invest 103869.61 1531586.9 221468.0 2045.8 value 1719503.68 27247303.7 3369944.3 38826.5 capital 413156.10 6153588.3 974281.3 8003.2 1 857.83 13418.2 1712.8 20.0 > crossprod(mmWH) invest value capital 1 invest 43732.40 643262.6 90592.41 857.83 value 643262.57 9942109.8 1344261.18 13418.20 capital 90592.41 1344261.2 220345.72 1712.80 1 857.83 13418.2 1712.80 20.00 > > ## p. 359 (higher precision than in Tab. 1) > coef(fmGE) (Intercept) value capital -9.95630645 0.02655119 0.15169387 > coef(fmWH) (Intercept) value capital -0.50939018 0.05289413 0.09240649 > > ## Scaled Sigma matrix for multivariate linear model (p. 359, middle), not reproducible > Y <- cbind(gr[1:20, 1], gr[21:40, 1]) > X <- as.matrix(cbind(gr[1:20, 3:2], 1, gr[21:40, 3:2], 1)) > B <- rbind( cbind(coef(fmGE)[3:1], 0), cbind(0, coef(fmWH)[3:1]) ) > crossprod((Y - X %*% B)) [,1] [,2] [1,] 13216.588 3528.981 [2,] 3528.981 1773.234 > > ## for comparison > fmSUR <- systemfit(invest ~ value + capital, data = pgr, method = "SUR", + control = systemfit.control(methodResidCov = "noDfCor")) > fmSUR$residCov * 20 General.Electric Westinghouse General.Electric 13788.376 3812.725 Westinghouse 3812.725 1801.301 > >