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.
> ## ----------------------------------------------------------------------------
> ## Koenker and Portnoy (1990)
> ##
> ## data
> ## availability: none
> ## firms: General Electric, Westinghouse
> ## errors: none
> ##
> ## analysis
> ## result: Table 1 - standard errors for SUR do not agree
> ## Table 2 - Westinghouse M estimation slightly different
> ## Table 3 - (NA)
> ## ----------------------------------------------------------------------------
>
> ## preliminaries
> 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
>
> ## data pre-processing
> gr <- subset(Grunfeld, firm %in% c("General Electric", "Westinghouse"))
> gr$value <- gr$value/1000
> gr$invest <- gr$invest/100
> gr$capital <- gr$capital/100
> gr$firm <- factor(gr$firm)
> pgr <- plm.data(gr, c("firm", "year"))
>
>
> ## Table 1, p. 1063
> ## OLS
> lm_OLS <- systemfit(invest ~ value + capital, method = "OLS", data = pgr)
> gsummary(lm_OLS)
Method: OLS
(Intercept) value capital R^2 sigma^2
General.Electric -0.100 0.266 0.152 0.705 0.078
0.314 0.156 0.026
Westinghouse -0.005 0.529 0.092 0.744 0.010
0.080 0.157 0.056
> ## alternatively:
> ## lm(invest ~ value + capital, data = gr, subset = firm == "General Electric")
> ## lm(invest ~ value + capital, data = gr, subset = firm == "Westinghouse")
>
> ## SUR
> ## these controls give the estimated "sigma matrix" and point estimates
> sf_SUR <- systemfit(invest ~ value + capital, method = "SUR", data = pgr,
+ control = systemfit.control(methodResidCov = "noDfCor"))
> gsummary(sf_SUR)
Method: SUR
(Intercept) value capital R^2 sigma^2
General.Electric -0.277 0.383 0.139 0.693 0.081
0.270 0.133 0.023
Westinghouse -0.013 0.576 0.064 0.740 0.011
0.070 0.134 0.049
> round(summary(sf_SUR)$residCov, digits = 3)
General.Electric Westinghouse
General.Electric 0.069 0.019
Westinghouse 0.019 0.009
>
> ## Table 2, p. 1064
> ## l_1 starting values
> library("quantreg")
Loading required package: SparseM
Package SparseM (0.83) loaded.
To cite, see citation("SparseM")
Attaching package: 'SparseM'
The following object(s) are masked from package:Matrix :
chol,
norm
The following object(s) are masked from package:stats :
model.response
The following object(s) are masked from package:base :
backsolve,
chol
Attaching package: 'quantreg'
The following object(s) are masked from package:survival :
untangle.specials
The following object(s) are masked from package:utils :
vignette
Warning message:
In getPackageName(env) :
Created a package name, "2010-01-29 18:21:43", when none found
> rq_GE <- rq(invest ~ value + capital, data = gr, subset = firm == "General Electric")
> rq_WH <- rq(invest ~ value + capital, data = gr, subset = firm == "Westinghouse")
> coef(rq_GE)
(Intercept) value capital
-0.1097989 0.2516002 0.1495661
> coef(rq_WH)
(Intercept) value capital
0.05076287 0.39702482 0.13927074
>
> ## M-estimation with custom psi function
> library("MASS")
> bread.rlm <- function(x, ...) {
+ xmat <- model.matrix(x)
+ xmat <- naresid(x$na.action, xmat)
+ wts <- weights(x)
+ if(is.null(wts)) wts <- 1
+ res <- residuals(x)
+ psi_deriv <- function(z) x$psi(z, deriv = 1)
+ rval <- sqrt(abs(as.vector(psi_deriv(res/x$s)/x$s))) * wts * xmat
+ rval <- chol2inv(qr.R(qr(rval))) * nrow(xmat)
+ return(rval)
+ }
>
> ## Koenker & Portnoy psi function
> psi.kp <- function(u, a = 0.99, b = 0.4, deriv = 0) {
+ lambda <- log((a + 1)/(1 - a)) / qnorm(1 - b)
+ if(!deriv) ifelse(u == 0, 0, (-(1 - 2/(1 + exp(-lambda * u))))/u)
+ else lambda * 2/(1 + exp(-lambda * u))^2 * exp(-lambda * u)
+ }
>
> ## Table 2, p. 1064
> ## M estimation (coefficients and standard errors)
> rlm_GE <- rlm(invest ~ value + capital, data = gr, subset = firm == "General Electric",
+ psi = psi.kp, init = coef(rq_GE), maxit = 100, acc = 1e-10)
> rlm_WH <- rlm(invest ~ value + capital, data = gr, subset = firm == "Westinghouse",
+ psi = psi.kp, init = coef(rq_WH), maxit = 100, acc = 1e-10)
> round(coeftest(rlm_GE, vcov = sandwich)[,1:2], digits = 3)
Estimate Std. Error
(Intercept) -0.119 0.072
value 0.252 0.028
capital 0.156 0.020
> round(coeftest(rlm_WH, vcov = sandwich)[,1:2], digits = 3)
Estimate Std. Error
(Intercept) 0.037 0.058
value 0.416 0.094
capital 0.135 0.041
>
> ## Table 3, p. 1064
> ## (would require some more coding)
>