knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
Load the hayashir
package
library(hayashir)
Let's get a quick look at our data by looking at the first 10 rows:
head(nerlove, 10)
Some information about the variables in the data can be found in the documentation:
?nerlove
And we can gain an understanding of the structure of our data by:
str(nerlove)
And by using the skim
command from the skimr
package we can look at summary statistics:
# install.packages("skimr") library(skimr) skim(nerlove)
unrestricted_ls <- log(total_cost) ~ log(output) + log(price_labor) + log(price_capital) + log(price_fuel) model1 <- lm(unrestricted_ls, data = nerlove) summary(model1)
Scale Coefficient:
return_to_scale <- 1 / model1$coefficients[2] print(return_to_scale)
t-stat:
The car
package provides the linearHypothesis
package which provides an easy way to test linear hypotheses:
library(car) linearHypothesis(model1, c("log(price_capital) = 0")) # relationship between t and F = F = t^2 test_results <- linearHypothesis(model1, c("log(price_capital) = 0")) # ie. same as above t_stat <- sqrt(test_results$F)
Following the book:
restricted_ls <- log(total_cost/price_fuel) ~ log(output) + log(price_labor/price_fuel) + log(price_capital / price_fuel) model2 <- lm(restricted_ls, data = nerlove) summary(model2)
Next, Hayashi provides a routine to compute the F-stat to test the restriction. First, we proceed as he instructs:
We can find SSR from the anova function:
anova(model1)
We need to get SSR_u from model 1 and the denominator df. The anova()
function returns a data.frame
from which we need to extract:
Sum_Sq
columnDf
columnTo get them:
anova_model1 <- anova(model1) final_row_u <- nrow(nrow(anova_model1)) SSR_u <- anova_model1$`Sum Sq`[final_row_u] df_resid_u <- anova_model1$Df[final_row_u]
and doing the same for the for the restricted model:
anova_model2 <- anova(model2) final_row_r <- nrow(nrow(anova_model2)) SSR_r <- anova_model2$`Sum Sq`[final_row_r] df_resid_r <- anova_model2$Df[final_row_r]
Then our F-stat can be computed as:
f_stat <- ( (SSR_r - SSR_u) / (df_resid_r - df_resid_u)) / (SSR_u / (df_resid_u)) print(paste("F stat is:", f_stat))
Which has a p-value of
pf(f_stat, df_resid_r - df_resid_u, df_resid_u, lower.tail = FALSE)
Alternatively, we can look at the F-stat versus a critical value.# The critical value at 5% is
qf(0.05, df_resid_r - df_resid_u, df_resid_u, lower.tail = FALSE)
Though instructive, that was kind of complicated ...
a simpler version would be using the linearHypothesis
function that we have already seen.
We specifying the restriction we want to impose:
linearHypothesis(model1, "log(price_labor) + log(price_capital) + log(price_fuel) = 1")
scale_effect <- log(total_cost / output) ~ log(output ) + log(price_labor ) + log(price_capital ) + log(price_fuel) model3 <- lm(scale_effect, data = nerlove) summary(model3)
linearHypothesis(model2, "log(output) = 1")
Which again returns an F-stat. We can map that into a t-stat as follows:
f_stat2 <- linearHypothesis(model2, "log(output) = 1") t_stat2 <- sqrt(f_stat2$F[2]) print(paste("t stat is:", t_stat2))
The car
package again helps us.
car
comes with a function residualPlot
which will plot residuals
against fitted values (by default),
or against a specified variable, in our case log(output)
# residual plot comes from car package residualPlot(model2, variable = "log(output)")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.