程序代写代做代考 econ61001_asn1.R

econ61001_asn1.R
2020-12-27

library(wooldridge)
## Warning: package ‘wooldridge’ was built under R version 4.0.3
library(MASS)
## 
## Attaching package: ‘MASS’
## The following object is masked from ‘package:wooldridge’:
## 
## cement
data(‘campus’)

##### Question 1
# calculate the sample mean of each explanatory variable, and round each value to 3 digits
round(apply(campus[, c(“enroll”, “police”, “priv”)], 2, mean), 3)
## enroll police priv 
## 16076.351 20.495 0.124
# calculate the sample standard deviation of each explanatory variable, and round each value to 3 digits
round(apply(campus[, c(“enroll”, “police”, “priv”)], 2, sd), 3)
## enroll police priv 
## 12298.986 15.631 0.331
##### Question 2
y <- campus$crime
x1 <- campus$enroll
x2 <- campus$police
# priv is a categorical variable
x3 <- as.factor(campus$priv)
# fitting the linear model using OLS
ols <- lm(y ~ x1+x2+x3)
summary(ols) ## 
## Call:
## lm(formula = y ~ x1 + x2 + x3)
## 
## Residuals:
## Min 1Q Median 3Q Max 
## -617.88 -123.67 1.68 82.41 1144.99 
## 
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) 
## (Intercept) -174.07758 47.67878 -3.651 0.000431 ***
## x1 0.02486 0.00290 8.572 2.14e-13 ***
## x2 7.78662 2.26744 3.434 0.000890 ***
## x31 75.05962 78.57454 0.955 0.341920 
## —
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ‘ 1
## 
## Residual standard error: 241.5 on 93 degrees of freedom
## Multiple R-squared: 0.7338, Adjusted R-squared: 0.7252 
## F-statistic: 85.46 on 3 and 93 DF, p-value: < 2.2e-16 # extract the residuals 
u <- resid(ols)
# x3^2 is omitted
f <- formula(u2 ~ (x1+x2+x3)^2+I(x1^2)+I(x2^2))
lm.fit <- lm(f, data=data.frame(u2=u^2, x1=x1, x2=x2, x3=x3))
summary(lm.fit) ## 
## Call:
## lm(formula = f, data = data.frame(u2 = u^2, x1 = x1, x2 = x2, 
## x3 = x3))
## 
## Residuals:
## Min 1Q Median 3Q Max 
## -230536 -55751 -10147 27014 1141347 
## 
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) 
## (Intercept) -3.867e+04 4.944e+04 -0.782 0.4362 
## x1 1.084e+01 4.983e+00 2.176 0.0323 *
## x2 -4.039e+03 4.351e+03 -0.928 0.3559 
## x31 4.177e+04 1.099e+05 0.380 0.7049 
## I(x1^2) -6.415e-06 1.150e-04 -0.056 0.9556 
## I(x2^2) 1.082e+02 7.753e+01 1.395 0.1664 
## x1:x2 -1.468e-01 1.464e-01 -1.003 0.3185 
## x1:x31 -6.174e+00 1.428e+01 -0.432 0.6665 
## x2:x31 1.872e+03 7.555e+03 0.248 0.8049 
## —
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ‘ 1
## 
## Residual standard error: 147000 on 88 degrees of freedom
## Multiple R-squared: 0.1893, Adjusted R-squared: 0.1156 
## F-statistic: 2.568 on 8 and 88 DF, p-value: 0.01446
n <- nrow(campus)
# the F-test statistic
F <- summary(lm.fit)$fstatistic
# the F-statistic
round(F[1], 3) ## value 
## 2.568 k <- F[2]
# the p-value for F-test version
round(pf(F[1],F[2],F[3], lower.tail = FALSE), 3) ## value 
## 0.014 # the R-squared of the regression
R2 <- summary(lm.fit)$r.squared
# the LM statistic
LM <- nrow(campus) * R2
# the LM statistic
round(LM, 3) ## [1] 18.358 # the p-value for LM-test version
round(pchisq(LM, df=k, lower.tail = FALSE), 3) ## [1] 0.019 ##### Question 3

mc <- 1000
F.result <- vector('logical', length = mc)
LM.result <- vector('logical', length = mc)
for (i in 1:mc) {
 set.seed(i)
 u <- rnorm(n)
 m <- lm(f,data=data.frame(u2=u^2, x1=x1, x2=x2, x3=x3))
 R2 <- summary(m)$r.squared
 F.result[i] <- (R2/k/((1-R2)/(n-k-1))) > F[1]
 LM.result[i] <- n*R2 > LM
}
# size F
sum(F.result)/mc
## [1] 0.05
# size LM
sum(LM.result)/mc
## [1] 0.05
Sigma <- diag(x1, n, n)
for (i in 1:mc) {
 set.seed(i)
 u <- mvrnorm(1, rep(0, n), Sigma)
 m <- lm(f,data=data.frame(u2=u^2, x1=x1, x2=x2, x3=x3))
 R2 <- summary(m)$r.squared
 F.result[i] <- (R2/k/((1-R2)/(n-k-1)))> F[1]
 LM.result[i] <- n*R2 > LM
}
# Power F D1
sum(F.result)/mc
## [1] 0.718
# Power LM D1
sum(LM.result)/mc
## [1] 0.718
Sigma <- diag(x1+x2, n, n)
for (i in 1:mc) {
 set.seed(i)
 u <- mvrnorm(1, rep(0, n), Sigma)
 m <- lm(f,data=data.frame(u2=u^2, x1=x1, x2=x2, x3=x3))
 R2 <- summary(m)$r.squared
 F.result[i] <- (R2/k/((1-R2)/(n-k-1)))> F[1]
 LM.result[i] <- n*R2 > LM
}
# Power F D2
sum(F.result)/mc
## [1] 0.713
# Power LM D2
sum(LM.result)/mc
## [1] 0.713