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)
>