knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Application: Returns to Scale in Electricity Supply

Load the hayashir package

library(hayashir)

The Data

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)

Testing Homogeneity of the cost function

The Unrestricted Model

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)

The Restricted Model

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:

  1. the SSR - located in the last row of the Sum_Sq column
  2. The degress of freedom - located in the last row of the Df column

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

Testing Restrictions in an Easier way

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

Detour: a cautionary note on $R^2$

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)

Testing Constant Returns to Scale

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

Importance of plotting residuals

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


lachlandeer/hayashir documentation built on Feb. 9, 2023, 2:01 p.m.