Empirical Project 7: Working in R code
Empirical Project 7 Working in R
These code downloads have been constructed as supplements to the full Doing Economics projects (https://core-econ.org/doing-economics/). You’ll need to download the data before running the code that follows.
Getting started in R
For this project you will need the following packages:
tidyverse, to help with data manipulationreadxl, to import an Excel spreadsheet.
If you need to install either of these packages, run the following code:
install.packages(c("readxl","tidyverse"))You can import the libraries now, or when they are used in the R walk-through below.
library(readxl)
library(tidyverse)## ── Attaching packages ──────────── tidyverse 1.3.0 ──## ✓ ggplot2 3.2.1 ✓ purrr 0.3.3
## ✓ tibble 2.1.3 ✓ dplyr 0.8.4
## ✓ tidyr 1.0.2 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.4.0## ── Conflicts ─────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()Part 7.1 Drawing supply and demand diagrams
R walk-through 7.1 Importing data into R and creating tables and charts
First we import the data with the read_excel function, using the na = "NA" option to indicate how missing data is recorded.
temp <- tempfile(fileext = ".xlsx")
dataURL <- "https://www.core-econ.org/index.php?uamfiletype=attachment&uamgetfile=/home/krzyszt2/public_html/wp-content/uploads/2018/06/Project-7-datafile.xlsx"
download.file(dataURL, destfile=temp, mode='wb')
wm_data <- readxl::read_excel(temp, sheet = "Sheet1", na = "NA")
str(wm_data)## Classes 'tbl_df', 'tbl' and 'data.frame': 22 obs. of 10 variables:
## $ Year : num 1930 1931 1932 1933 1934 ...
## $ log q
## (Q) : num 4.45 4.36 4.2 4.03 4.1 ...
## $ log h
## (X) : num 4.38 4.33 4.05 4.01 4.09 ...
## $ log p
## (P) : num 4.76 4.61 4.37 4.53 4.64 ...
## $ log p_c
## (C) : num 2.25 1.73 1.87 2.32 2.51 ...
## $ log p_v
## (T) : num 0.845 2.726 2.588 2.286 1.476 ...
## $ log w
## (W) : num 3.37 3.14 2.83 2.77 2.92 ...
## $ log n
## (N) : num 4.81 4.82 4.83 4.83 4.84 ...
## $ log(y/n)
## (Y/N): num 6.4 6.24 5.97 5.9 6.02 ...
## $ log p_f
## (F) : num 2.54 2.55 2.6 2.65 2.62 ...Let’s use the exp function to create the variables p and q from their log counterparts (renamed as log.p and log.q respectively). We also transform the harvest variable (renamed as log.h) and save it as h. The harvest will be at most as large as the crop (q).
names(wm_data)## [1] "Year" "log q \r\n(Q)" "log h \r\n(X)"
## [4] "log p \r\n(P)" "log p_c \r\n(C)" "log p_v \r\n(T)"
## [7] "log w \r\n(W)" "log n \r\n(N)" "log(y/n) \r\n(Y/N)"
## [10] "log p_f \r\n(F)"names(wm_data) <- c("Year", "log.q", "log.h", "log.p", "log.pc", "log.pv", "log.w", "log.n", "log.yn", "log.pf")
wm_data <- wm_data %>% mutate(p = exp(log.p), # Price
h = exp(log.h), # Harvest quantity
q = exp(log.q)) # Crop quantityLet’s use ggplot to produce the chart for the prices, with Year as the horizontal axis variable (xlab) and price (p) as the vertical axis variable (ylab).
wm_data %>% ggplot(aes(x=Year, y=p)) + geom_line() + geom_point() +
labs(title = "Price of Watermelon") + xlab("Year") + ylab("Price")
Now we create the line chart for harvest and crop quantities (the variables h and q, respectively). First, we plot the crop quantities as a dashed line (lty = "dashed"), then use lines to add a solid line for the harvest data. The legend function adds a chart legend at the specified coordinates (the first two arguments in the function).
wm_data %>% select(Year, q, h) %>%
rename(Crop = q,
Harvest = h) %>%
pivot_longer(-Year, names_to = "Activity") %>%
ggplot(aes(x = Year, y = value, group = Activity)) +
geom_line(aes(linetype = Activity, color = Activity)) +
geom_point(aes(color = Activity)) +
labs(title = "Crop and Harvest of Watermelon") + xlab("Year") + ylab("Millions")
[End of walk-through]
log(P) = -2 + 1.7 log(Q) (Supply curve)
log(P) = 8.5 - 0.82 log(Q) (Demand curve)
table <- tibble(Q=seq(20,100,5))
table <- table %>% mutate(lnQ=log(Q),
logP_supply=-2+1.7*lnQ,
logP_demand=8.5-0.82*lnQ,
Supply=exp(logP_supply),
Demand=exp(logP_demand))
tablePlot of supply and demand functions.
table %>% select(Q, Demand, Supply) %>%
pivot_longer(-Q, names_to = "Curve", values_to = "P") %>%
ggplot(aes(x=Q, y=P, color=Curve)) + geom_line() + ggtitle("Supply and Demand for Watermelons")
Plot of new supply shift and the demand function.
table <- table %>% mutate(New_Supply=exp(-2+1.7*lnQ+0.4))
table %>% select(Q, Demand, Supply, New_Supply) %>%
pivot_longer(-Q, names_to = "Curve", values_to = "P") %>%
ggplot(aes(x=Q, y=P, color=Curve)) + geom_line() +
labs(title = "Supply and Demand for Watermelons",
subtitle = "with supply shift")
Surpluses
We can now work with the functions directly to calculate the consumer and producer surplus. We start by finding the equilibrium quantity (\(q_1\)). The equilibrium quantity (“x” in the formulas) is found by setting the supply function equal to demand function, and then solve for quantity (\(q_1\)), which is at the intersection of the supply and demand curve.
suppressPackageStartupMessages(library(mosaic))
suppressPackageStartupMessages(library(mosaicCalc))
supply <- makeFun(exp(-2+1.7*log(x))~x)
demand <- makeFun(exp(8.5-0.82*log(x))~x)
q1 <- findZeros(supply(x)-demand(x) ~ x, x.lim=range(60,80))
q1 <- q1$x
q1## [1] 64.5001We then insert the equilibrium quantity in the supply curve (could also have used the demand curve) to find the equilibrium price (\(p_1\)).
p1 <- supply(q1)
p1## [1] 161.3109Let’s update the plot with the equilibrium values in it:
Q <- 20:100
ggplot() +
stat_function(aes(x=Q, color = "Demand"), fun = demand) +
stat_function(aes(x=Q, color = "Supply"), fun = supply) +
labs(x="Quantity", y="Price", color=NULL) +
annotate("point", x = q1, y = p1, color = "grey30") +
annotate("segment", x = q1, xend = q1, y = 20, yend = p1,
linetype = "dashed", color = "grey30") +
annotate("segment", x = 20, xend = q1, y = p1, yend = p1,
linetype = "dashed", color = "grey30") +
annotate("text", x = 66, y = 10, label = q1, parse = TRUE) +
annotate("text", x = 25, y = 170, label = p1, parse = TRUE)
The numerical values of the surpluses are found by the following relationships:
consumer surplus = \(\int_{q_0}^{q_1}D(x)dx - p_1 \times (q_1-q_0)\)
where \(q_0\) is the starting value of quantity, in our case 20, and \(q_1\) is the equilibrium quantity.
Similarly, we find the producer surlus as:
producer surplus = \(p_1 \times (q_1-q_0) - \int_{q_0}^{q_1}S(x)dx\)
Before we do the calculations, we need to find the anti derivative of the supply and demand function in order to be able to calculate the itegrals.
S <- antiD(supply(x)~x)
D <- antiD(demand(x)~x)Now, consumer surplus calculated in the quantity interval {20, to equilibrium quantity: 64.5001} is:
(D(q1)-D(20))-p1*(q1-20)## [1] 3806.455Producer surplus calculated in the quantity interval {20, to equilibrium quantity: 64.5001} is:
p1*(q1-20)-(S(q1)-S(20))## [1] 3488.048The total surplus is consumer surplus plus producer surplus:
(D(q1)-D(20))-p1*(q1-20)+p1*(q1-20)-(S(q1)-S(20))## [1] 7294.503Another way of calculating the total surplus is to take the integral under the demand curve, and subtract the integral under the supply curve:
(D(q1)-D(20))-(S(q1)-S(20))## [1] 7294.503Let’s make a plot of these areas, and format the plot as used in “The Economist”.
z <- seq(20, q1, 0.01)
ggplot() +
stat_function(aes(x=Q, color = "Demand"), fun = demand) +
stat_function(aes(x=Q, color = "Supply"), fun = supply) +
labs(x="Quantity", y="Price", color=NULL, fill=NULL) +
annotate("point", x = q1, y = p1, color = "grey30") +
annotate("segment", x = q1, xend = q1, y = 20, yend = p1,
linetype = "dashed", color = "grey30") +
annotate("segment", x = 20, xend = q1, y = p1, yend = p1,
linetype = "dashed", color = "grey30") +
geom_ribbon(aes(x = z, ymin = supply(z), ymax = p1,
fill = "Producer surplus"), alpha = 0.25) +
geom_ribbon(aes(x = z, ymin = p1, ymax = demand(z),
fill = "Consumer surplus"), alpha = 0.25) + ggthemes::theme_economist()
The following function utilizes the predefined supply and demand functions as inputs to find the equilibrium quantity and price, and then calculates the surpluses. One also has to specify the quantity domain in which the calculations take place, the smallest value is \(q_0\), while the larger value has to be \(\geqslant\) than the equilibrium quantity \(q_1\), since the function calculates the equilibrium quantity (and price) on the fly.
surpluses <- function(demand, supply, qdomain) {
q <- uniroot(function(x) demand(x) - supply(x), qdomain)$root
p <- supply(q)
consumer_surplus <- integrate(demand, qdomain[1], q)$value - p*(q-qdomain[1])
producer_surplus <- p*(q-qdomain[1]) - integrate(supply, qdomain[1], q)$value
total_surplus <- consumer_surplus + producer_surplus
list("consumer" = consumer_surplus,
"producer" = producer_surplus,
"total" = total_surplus,
"eq_p" = p,
"eq_q" = q) }
surpluses(demand, supply, c(20, 100))## $consumer
## [1] 3806.459
##
## $producer
## [1] 3488.044
##
## $total
## [1] 7294.503
##
## $eq_p
## [1] 161.3108
##
## $eq_q
## [1] 64.50008New Supply
The new supply curve is:
new_supply <- makeFun(exp(-2+1.7*log(x)+0.4)~x)With the following surpluses:
surpluses(demand, new_supply, c(20, 100))## $consumer
## [1] 2919.849
##
## $producer
## [1] 2935.335
##
## $total
## [1] 5855.184
##
## $eq_p
## [1] 183.7346
##
## $eq_q
## [1] 55.03319This indicates a reduction in total surplus of:
old <- surpluses(demand, supply, c(20, 100))
new <- surpluses(demand, new_supply, c(20, 100))
old$total-new$total## [1] 1439.319Producer surplus is reduced by:
old$producer-new$producer## [1] 552.7093or 15.85 percent.
Consumer surplus is reduced by:
old$consumer-new$consumer## [1] 886.6101or 23.29 percent.
Since the demand curve is “more” inelastic than the supply curve, more of the reduction in surplus is on the consumer side.
Estimating supply and demand
names(wm_data)## [1] "Year" "log.q" "log.h" "log.p" "log.pc" "log.pv" "log.w" "log.n"
## [9] "log.yn" "log.pf" "p" "h" "q"wm_data <- wm_data %>% mutate(CP = ifelse(Year %in% 1934:1951,1,0),
WW2 = ifelse(Year %in% 1943:1946,1,0))
supply_fit <- lm(log.q ~ dplyr::lag(log.p) + dplyr::lag(log.pc) + dplyr::lag(log.pv) + CP + WW2, data = wm_data)
summary(supply_fit)##
## Call:
## lm(formula = log.q ~ dplyr::lag(log.p) + dplyr::lag(log.pc) +
## dplyr::lag(log.pv) + CP + WW2, data = wm_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0039962 -0.0015029 0.0001659 0.0013821 0.0046611
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.4240906 0.0088665 273.40 < 2e-16 ***
## dplyr::lag(log.p) 0.5796588 0.0031290 185.26 < 2e-16 ***
## dplyr::lag(log.pc) -0.3219423 0.0030646 -105.05 < 2e-16 ***
## dplyr::lag(log.pv) -0.1235780 0.0007502 -164.73 < 2e-16 ***
## CP 0.0729619 0.0019863 36.73 4.15e-16 ***
## WW2 -0.3601240 0.0021652 -166.32 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.002473 on 15 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.9998, Adjusted R-squared: 0.9997
## F-statistic: 1.375e+04 on 5 and 15 DF, p-value: < 2.2e-16round(coef(supply_fit),3)## (Intercept) dplyr::lag(log.p) dplyr::lag(log.pc) dplyr::lag(log.pv)
## 2.424 0.580 -0.322 -0.124
## CP WW2
## 0.073 -0.360round(confint(supply_fit),3)## 2.5 % 97.5 %
## (Intercept) 2.405 2.443
## dplyr::lag(log.p) 0.573 0.586
## dplyr::lag(log.pc) -0.328 -0.315
## dplyr::lag(log.pv) -0.125 -0.122
## CP 0.069 0.077
## WW2 -0.365 -0.356Slight deviation between dataset values and values in Wold. (https://www.cebm.net/2019/03/logarithms-and-log-transformations-2/)
# last value of log q column of data in natural logs
log(10^1.921)## [1] 4.423266Actual value in log10
wm_data$log.q[22]*log10(exp(1))## [1] 1.921012should be
10^1.921## [1] 83.36812it is
10^(wm_data$log.q[22]*log10(exp(1)))## [1] 83.37051Fix:
as.numeric(substr(wm_data$log.q*log10(exp(1)),1,5))## [1] 1.932 1.892 1.826 1.751 1.779 1.822 1.796 1.851 1.850 1.800 1.887 1.809
## [13] 1.752 1.686 1.850 1.863 1.908 1.909 1.868 1.894 1.916 1.921Estimate the demand model using OLS.
demand_fit <- lm(log.p ~ I(log.h-log.n) + log.yn + log.pf, data = wm_data)
summary(demand_fit)##
## Call:
## lm(formula = log.p ~ I(log.h - log.n) + log.yn + log.pf, data = wm_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.28656 -0.07373 0.01239 0.08581 0.18532
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.63865 0.46632 -7.803 3.49e-07 ***
## I(log.h - log.n) -0.94372 0.22780 -4.143 0.000611 ***
## log.yn 1.55449 0.08712 17.843 6.84e-13 ***
## log.pf -0.73953 0.18617 -3.972 0.000893 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1308 on 18 degrees of freedom
## Multiple R-squared: 0.9626, Adjusted R-squared: 0.9564
## F-statistic: 154.4 on 3 and 18 DF, p-value: 4.957e-13logp <- as.numeric(substr(wm_data$log.p*log10(exp(1)),1,5))
logh <- as.numeric(substr(wm_data$log.h*log10(exp(1)),1,5))
logn <- as.numeric(substr(wm_data$log.n*log10(exp(1)),1,5))
logyn <- as.numeric(substr(wm_data$log.yn*log10(exp(1)),1,5))
logpf <- as.numeric(substr(wm_data$log.pf*log10(exp(1)),1,5))
summary(lm(logp ~ I(logh-logn) + logyn + logpf))##
## Call:
## lm(formula = logp ~ I(logh - logn) + logyn + logpf)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.124449 -0.032019 0.005382 0.037266 0.080482
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.58024 0.20252 -7.803 3.49e-07 ***
## I(logh - logn) -0.94372 0.22780 -4.143 0.000611 ***
## logyn 1.55449 0.08712 17.843 6.84e-13 ***
## logpf -0.73953 0.18617 -3.972 0.000893 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0568 on 18 degrees of freedom
## Multiple R-squared: 0.9626, Adjusted R-squared: 0.9564
## F-statistic: 154.4 on 3 and 18 DF, p-value: 4.957e-13The intercept changes.
wm_data <- wm_data %>% mutate(Llnp = dplyr::lag(log.p),
Llnpc = dplyr::lag(log.pc),
Llnpv = dplyr::lag(log.pv),
log.hn = log.h-log.n,
log.pw = log.p-log.w,
log.y = log.yn + log.n)
names(wm_data)## [1] "Year" "log.q" "log.h" "log.p" "log.pc" "log.pv" "log.w" "log.n"
## [9] "log.yn" "log.pf" "p" "h" "q" "CP" "WW2" "Llnp"
## [17] "Llnpc" "Llnpv" "log.hn" "log.pw" "log.y"library(systemfit)## Loading required package: car## Loading required package: carData##
## Attaching package: 'car'## The following object is masked _by_ '.GlobalEnv':
##
## S## The following object is masked from 'package:mosaicCore':
##
## logit## The following objects are masked from 'package:mosaic':
##
## deltaMethod, logit## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:purrr':
##
## some## Loading required package: lmtest## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric##
## Please cite the 'systemfit' package as:
## Arne Henningsen and Jeff D. Hamann (2007). systemfit: A Package for Estimating Systems of Simultaneous Equations in R. Journal of Statistical Software 23(4), 1-40. http://www.jstatsoft.org/v23/i04/.
##
## If you have questions, suggestions, or comments regarding the 'systemfit' package, please use a forum or 'tracker' at systemfit's R-Forge site:
## https://r-forge.r-project.org/projects/systemfit/# Specify the three equations
cropSupply <- log.q ~ Llnp + Llnpc + Llnpv + CP + WW2
harvestSupply <- log.h ~ log.pw + log.q
demand <- log.p ~ log.hn + log.yn + log.pf
# Put them in a list
eqSystem <- list(demand = demand, supply1 = cropSupply, supply2 = harvestSupply)
# Estimate the system
fit2sls <- systemfit(eqSystem, method = "2SLS", inst=~Llnp+log.w+log.n+log.y+log.pf+CP+WW2+Llnpc+Llnpv, data = wm_data)
print(fit2sls)##
## systemfit results
## method: 2SLS
##
## Coefficients:
## demand_(Intercept) demand_log.hn demand_log.yn demand_log.pf
## -3.4494342 -0.9085628 1.5634243 -0.8194997
## supply1_(Intercept) supply1_Llnp supply1_Llnpc supply1_Llnpv
## 2.4246040 0.5792040 -0.3209264 -0.1237789
## supply1_CP supply1_WW2 supply2_(Intercept) supply2_log.pw
## 0.0726426 -0.3601231 -0.5232387 0.1254032
## supply2_log.q
## 1.0673956summary(fit2sls)##
## systemfit results
## method: 2SLS
##
## N DF SSR detRCov OLS-R2 McElroy-R2
## system 60 47 0.316663 0 0.961725 0.999282
##
## N DF SSR MSE RMSE R2 Adj R2
## demand 20 16 0.294269 0.018392 0.135616 0.960975 0.953657
## supply1 20 14 0.000088 0.000006 0.002508 0.999770 0.999688
## supply2 20 17 0.022306 0.001312 0.036223 0.936170 0.928661
##
## The covariance matrix of the residuals
## demand supply1 supply2
## demand 0.018391786 7.39640e-05 5.46955e-04
## supply1 0.000073964 6.28796e-06 -5.72385e-06
## supply2 0.000546955 -5.72385e-06 1.31213e-03
##
## The correlations of the residuals
## demand supply1 supply2
## demand 1.000000 0.2174972 0.1113400
## supply1 0.217497 1.0000000 -0.0630152
## supply2 0.111340 -0.0630152 1.0000000
##
##
## 2SLS estimates for 'demand' (equation 1)
## Model Formula: log.p ~ log.hn + log.yn + log.pf
## Instruments: ~Llnp + log.w + log.n + log.y + log.pf + CP + WW2 + Llnpc + Llnpv
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.4494342 0.5649728 -6.10549 1.5186e-05 ***
## log.hn -0.9085628 0.2776009 -3.27291 0.0047846 **
## log.yn 1.5634243 0.0909038 17.19867 9.6660e-12 ***
## log.pf -0.8194997 0.2215706 -3.69859 0.0019481 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.135616 on 16 degrees of freedom
## Number of observations: 20 Degrees of Freedom: 16
## SSR: 0.294269 MSE: 0.018392 Root MSE: 0.135616
## Multiple R-Squared: 0.960975 Adjusted R-Squared: 0.953657
##
##
## 2SLS estimates for 'supply1' (equation 2)
## Model Formula: log.q ~ Llnp + Llnpc + Llnpv + CP + WW2
## Instruments: ~Llnp + log.w + log.n + log.y + log.pf + CP + WW2 + Llnpc + Llnpv
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.424603997 0.009017399 268.8806 < 2.22e-16 ***
## Llnp 0.579203972 0.003228746 179.3898 < 2.22e-16 ***
## Llnpc -0.320926429 0.003380762 -94.9273 < 2.22e-16 ***
## Llnpv -0.123778900 0.000805006 -153.7615 < 2.22e-16 ***
## CP 0.072642607 0.002057384 35.3082 4.4409e-15 ***
## WW2 -0.360123062 0.002195940 -163.9949 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.002508 on 14 degrees of freedom
## Number of observations: 20 Degrees of Freedom: 14
## SSR: 8.8e-05 MSE: 6e-06 Root MSE: 0.002508
## Multiple R-Squared: 0.99977 Adjusted R-Squared: 0.999688
##
##
## 2SLS estimates for 'supply2' (equation 3)
## Model Formula: log.h ~ log.pw + log.q
## Instruments: ~Llnp + log.w + log.n + log.y + log.pf + CP + WW2 + Llnpc + Llnpv
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.5232387 0.4800285 -1.09002 0.29092
## log.pw 0.1254032 0.0611267 2.05153 0.05595 .
## log.q 1.0673956 0.0936429 11.39858 2.2005e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.036223 on 17 degrees of freedom
## Number of observations: 20 Degrees of Freedom: 17
## SSR: 0.022306 MSE: 0.001312 Root MSE: 0.036223
## Multiple R-Squared: 0.93617 Adjusted R-Squared: 0.928661# Estimated as a simultaneous system, like 2SLS + SUR
fit3sls <- systemfit(eqSystem, method = "3SLS", inst=~Llnp+log.w+log.n+log.y+log.pf+CP+WW2+Llnpc+Llnpv, data = wm_data)
summary(fit3sls)##
## systemfit results
## method: 3SLS
##
## N DF SSR detRCov OLS-R2 McElroy-R2
## system 60 47 0.31732 0 0.961646 0.999297
##
## N DF SSR MSE RMSE R2 Adj R2
## demand 20 16 0.294973 0.018436 0.135779 0.960881 0.953547
## supply1 20 14 0.000089 0.000006 0.002521 0.999768 0.999685
## supply2 20 17 0.022258 0.001309 0.036184 0.936307 0.928814
##
## The covariance matrix of the residuals used for estimation
## demand supply1 supply2
## demand 0.018391786 7.39640e-05 5.46955e-04
## supply1 0.000073964 6.28796e-06 -5.72385e-06
## supply2 0.000546955 -5.72385e-06 1.31213e-03
##
## The covariance matrix of the residuals
## demand supply1 supply2
## demand 1.84358e-02 9.09893e-05 6.15266e-04
## supply1 9.09893e-05 6.35748e-06 -7.27603e-06
## supply2 6.15266e-04 -7.27603e-06 1.30930e-03
##
## The correlations of the residuals
## demand supply1 supply2
## demand 1.000000 0.2657766 0.1252311
## supply1 0.265777 1.0000000 -0.0797502
## supply2 0.125231 -0.0797502 1.0000000
##
##
## 3SLS estimates for 'demand' (equation 1)
## Model Formula: log.p ~ log.hn + log.yn + log.pf
## Instruments: ~Llnp + log.w + log.n + log.y + log.pf + CP + WW2 + Llnpc + Llnpv
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.3159607 0.5591173 -5.93071 2.1121e-05 ***
## log.hn -0.8935079 0.2766901 -3.22927 0.0052448 **
## log.yn 1.5590872 0.0903938 17.24773 9.2559e-12 ***
## log.pf -0.8536652 0.2190359 -3.89738 0.0012808 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.135779 on 16 degrees of freedom
## Number of observations: 20 Degrees of Freedom: 16
## SSR: 0.294973 MSE: 0.018436 Root MSE: 0.135779
## Multiple R-Squared: 0.960881 Adjusted R-Squared: 0.953547
##
##
## 3SLS estimates for 'supply1' (equation 2)
## Model Formula: log.q ~ Llnp + Llnpc + Llnpv + CP + WW2
## Instruments: ~Llnp + log.w + log.n + log.y + log.pf + CP + WW2 + Llnpc + Llnpv
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.424497799 0.008910393 272.0977 < 2.22e-16 ***
## Llnp 0.579253472 0.003178878 182.2194 < 2.22e-16 ***
## Llnpc -0.321131860 0.003303094 -97.2215 < 2.22e-16 ***
## Llnpv -0.123731712 0.000797715 -155.1077 < 2.22e-16 ***
## CP 0.073081596 0.002001878 36.5065 2.6645e-15 ***
## WW2 -0.360617909 0.002177444 -165.6153 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.002521 on 14 degrees of freedom
## Number of observations: 20 Degrees of Freedom: 14
## SSR: 8.9e-05 MSE: 6e-06 Root MSE: 0.002521
## Multiple R-Squared: 0.999768 Adjusted R-Squared: 0.999685
##
##
## 3SLS estimates for 'supply2' (equation 3)
## Model Formula: log.h ~ log.pw + log.q
## Instruments: ~Llnp + log.w + log.n + log.y + log.pf + CP + WW2 + Llnpc + Llnpv
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.4956091 0.4796622 -1.03325 0.315965
## log.pw 0.1215581 0.0610608 1.99077 0.062831 .
## log.q 1.0623814 0.0935825 11.35235 2.3401e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.036184 on 17 degrees of freedom
## Number of observations: 20 Degrees of Freedom: 17
## SSR: 0.022258 MSE: 0.001309 Root MSE: 0.036184
## Multiple R-Squared: 0.936307 Adjusted R-Squared: 0.928814print(fit2sls)##
## systemfit results
## method: 2SLS
##
## Coefficients:
## demand_(Intercept) demand_log.hn demand_log.yn demand_log.pf
## -3.4494342 -0.9085628 1.5634243 -0.8194997
## supply1_(Intercept) supply1_Llnp supply1_Llnpc supply1_Llnpv
## 2.4246040 0.5792040 -0.3209264 -0.1237789
## supply1_CP supply1_WW2 supply2_(Intercept) supply2_log.pw
## 0.0726426 -0.3601231 -0.5232387 0.1254032
## supply2_log.q
## 1.0673956print(fit3sls)##
## systemfit results
## method: 3SLS
##
## Coefficients:
## demand_(Intercept) demand_log.hn demand_log.yn demand_log.pf
## -3.3159607 -0.8935079 1.5590872 -0.8536652
## supply1_(Intercept) supply1_Llnp supply1_Llnpc supply1_Llnpv
## 2.4244978 0.5792535 -0.3211319 -0.1237317
## supply1_CP supply1_WW2 supply2_(Intercept) supply2_log.pw
## 0.0730816 -0.3606179 -0.4956091 0.1215581
## supply2_log.q
## 1.0623814length(coef(fit3sls))## [1] 13library(car)
deltaMethod(fit3sls, "1/demand_log.hn") # own-price elasticitydeltaMethod(fit3sls, "-demand_log.yn/demand_log.hn")deltaMethod(fit3sls, "-demand_log.pf/demand_log.hn")R LIML estimation (https://cran.r-project.org/web/packages/ivmodel/ivmodel.pdf)