# 34.76626 47.59049
y_min <- min(sortBy_pA$Q_dot_pA)
y_max <- max(sortBy_pA$Q_dot_pA)
print(cbind(y_min, y_max))
# 20548.31 26708.34
print(sortBy_pA)
P_apple_asX <- data.frame(P_apple = sortBy_pA$P_apple)
print(P_apple_asX)
Q_mix_asY <- data.frame(Q_fit_aA = sortBy_pA$Q_fit_pA, 
                        Q_dot_pA = sortBy_pA$Q_dot_pA)
print(Q_mix_asY)
plot.new()
matplot(P_apple_asX, 
        Q_mix_asY, 
  type = c("l", "p"), 
  lwd = 2, pch = 20, col = 1, 
  xlim = c(34.0, 48.0), 
  ylim = c(20000, 28000), 
  xlab = "apple price [yen/100g]", 
  ylab = "apple quantity [g]", 
  main = "inverse per-capita demand curve in raw scale")
#### additional infromation
# text(x_max, Q_mix_asY$Q_fit_aA[P_apple_asX == x_max], 
#    labels = "invD_apple", pos = 4)
#### save the graph as a file
# dev.copy(jpeg, file = "outputs/inv_demand_curve_rawScale_wDot.jpg")
# dev.off()
#### plotting P as y-axis, Q as x-axis with plot()
#### sort by Q_apple
sortBy_qA <- xy_mix[order(xy_mix$Q_fit_pA), ]
Q_fit_asX <- sortBy_qA$Q_fit_pA
P_apple_asY <- sortBy_qA$P_apple
x_min <- min(Q_fit_asX)
x_max <- max(Q_fit_asX)
print(cbind(x_min, x_max))
# 20830.26 25118.47
y_max <- max(P_apple_asY)
y_min <- min(P_apple_asY)
print(cbind(y_min, y_max))
# 34.76626 47.59049
#### demand curve in raw scale
plot.new()
plot(Q_fit_asX, P_apple_asY, 
  type = "l", lwd = 2, col = 1, 
  xlim = c(20000, 27000), 
  ylim = c(34.0, 48.0), 
  xlab = "apple quantity [g]", 
  ylab = "apple price [yen/100g]", 
  main = "demand curve in raw scale")
#### additional information
text(x_max, P_apple_asY[Q_fit_asX == x_max], 
    labels = "D_apple", pos = 4)
#### save the graph into a file
# dev.copy(jpeg, file = "./outputs/graph_demand_curve_apple.jpg")
# dev.off()
### end #############################################
### R Script ########################################
# Last modfied: 2025-07-03
# Name of script: graph_demand_curve_apple.R
# Keywords: 
# graph_demand_curve: get the graph of demand curve,
# apple: the demand curve of apples,
#
# Written by: Hiroshi Murao
#
# Purpose: regression analysis about apple demand in Japan. 
#
# Original data source: ƌvNñ published by ӝvÁE
# The data can be obtained by visiting the following Web sites.
# ӝvǁAĉƌvN??Eҍ،ˉʁAēѐlȏÂ̐сAĕݐ̂PONxN??E݌v\ADƂPV\@
# H׿u̍؁ECєiжN̍؁jv`uʕȁv
#
# The following is a list of variable names in the data set. 
#
# year
# EXP_fresh_fruitQ_fresh_fruitPN_fresh_fruitP_fresh_fruit
# EXP_appleQ_applePN_appleP_apple
# EXP_mandarinQ_mandarinPN_mandarinP_mandarin
# EXP_pearQ_pearPN_pearP_pear
# EXP_grapeQ_grapePN_grapeP_grape
# EXP_persimmonQ_pwersimmonPN_persimmonP_persimmon
# EXP_peachQ_peachPN_peachP_peach
# EXP_watermelonQ_watermelonPN_watermelonP_watermelon
# EXP_melonQ_melonPN_melonP_melon
# EXP_strawberryQ_strawberryPN_strawberryP_strawberry
# EXP_bananaQ_bananaPN_bananaP_banana
# EXP_other_fruitQ_other_fruitPN_other_fruitP_other_fruit
# CPI
# income
# trend
#
# Input: data of the above variables.
# Output: estimation results and estimated curves
#
##########################################################
### first things to do 
rm(list = ls())
setwd("C:/book_SIMU/data_works/Apple/R_prg")
### read data
in_data <- read.csv("./data/data_apple_JP.csv")
names(in_data)
# print(in_data)
### estimate the apple demand curve
out_lm = lm(log(Q_apple) ~ log(P_apple) + log(income) + log(P_persimmon) 
          + trend, data = in_data)
trend2pw <- in_data$trend^2
# out_lm = lm(log(Q_apple) ~ log(P_apple) + log(income) + log(P_persimmon) 
#            + trend + trend2pw, data = in_data)   
summary(out_lm)
#### preparing for plotting curves
b <- coefficients(out_lm)
resid <- residuals(out_lm)
# print(b)
logQ_fit_pA <- b[1] + 
               b[2]*(log(in_data$P_apple)) + 
               b[3]*mean(log(in_data$income)) +  
               b[4]*mean(log(in_data$P_persimmon) + 
               b[5]*mean(in_data$trend))
Q_fit_pA <- exp(logQ_fit_pA)
logQ_dot_pA <- logQ_fit_pA + resid
Q_dot_pA <- exp(logQ_dot_pA)
#### plot: logQ as y-axis, logP as x-axis, with matplot()
logP_apple <- log(in_data$P_apple)
print(logP_apple)
x_min <- min(logP_apple) 
x_max <- max(logP_apple)
print(cbind(x_min, x_max))
# 3.548647 3.862633
logQ_mix <- data.frame(logQ_fit_pA = logQ_fit_pA, 
                       logQ_dot_pA = logQ_dot_pA)
y_min <- min(logQ_mix$logQ_dot_pA) 
y_max <- max(logQ_mix$logQ_dot_pA)
print(cbind(y_min, y_max))
# 9.930534 10.19273
#### inverse demand curve
plot.new()
matplot(logP_apple, logQ_mix, 
  type = c("l", "p"), 
  lwd = 2, pch = 20, col = 1, 
  xlim = c(3.5, 3.9), 
  ylim = c(9.9, 10.2), 
  xlab = "log of apple price [yen/100g]", 
  ylab = "log of apple quantity [g]", 
  main = "inverse demand curve in log scale")
#### additional infromation
# text(x_max, logQ_mix$logQ_fit_aA[logP_apple == x_max], 
#    labels = "invD_apple", pos = 4)
#### save as a file
# dev.copy(jpeg, file = "./outputs/inv_demand_curve_logScale.jpg")
# dev.off()
#### plotting Q as y-axis, P as x-axis, with matplot()
xy_mix <- data.frame(P_apple = in_data$P_apple, 
                     Q_fit_pA = Q_fit_pA, 
                     Q_dot_pA = Q_dot_pA) 
print(xy_mix)
#### sort by P_apple
sortBy_pA <- xy_mix[order(xy_mix$P_apple), ]
x_min <- min(sortBy_pA$P_apple)
x_max <- max(sortBy_pA$P_apple)
print(cbind(x_min, x_max))
# 34.76626 47.59049
y_min <- min(sortBy_pA$Q_dot_pA)
y_max <- max(sortBy_pA$Q_dot_pA)
print(cbind(y_min, y_max))
# 20548.31 26708.34
print(sortBy_pA)
P_apple_asX <- data.frame(P_apple = sortBy_pA$P_apple)
print(P_apple_asX)
Q_mix_asY <- data.frame(Q_fit_aA = sortBy_pA$Q_fit_pA, 
                        Q_dot_pA = sortBy_pA$Q_dot_pA)
print(Q_mix_asY)
plot.new()
matplot(P_apple_asX, 
        Q_mix_asY, 
  type = c("l", "p"), 
  lwd = 2, pch = 20, col = 1, 
  xlim = c(34.0, 48.0), 
  ylim = c(20000, 28000), 
  xlab = "apple price [yen/100g]", 
  ylab = "apple quantity [g]", 
  main = "inverse per-capita demand curve in raw scale")
#### additional infromation
# text(x_max, Q_mix_asY$Q_fit_aA[P_apple_asX == x_max], 
#    labels = "invD_apple", pos = 4)
#### save the graph as a file
# dev.copy(jpeg, file = "outputs/inv_demand_curve_rawScale_wDot.jpg")
# dev.off()
#### plotting P as y-axis, Q as x-axis with plot()
#### sort by Q_apple
sortBy_qA <- xy_mix[order(xy_mix$Q_fit_pA), ]
Q_fit_asX <- sortBy_qA$Q_fit_pA
P_apple_asY <- sortBy_qA$P_apple
x_min <- min(Q_fit_asX)
x_max <- max(Q_fit_asX)
print(cbind(x_min, x_max))
# 20830.26 25118.47
y_max <- max(P_apple_asY)
y_min <- min(P_apple_asY)
print(cbind(y_min, y_max))
# 34.76626 47.59049
#### demand curve in raw scale
plot.new()
plot(Q_fit_asX, P_apple_asY, 
  type = "l", lwd = 2, col = 1, 
  xlim = c(20000, 27000), 
  ylim = c(34.0, 48.0), 
  xlab = "apple quantity [g]", 
  ylab = "apple price [yen/100g]", 
  main = "demand curve in raw scale")
#### additional information
text(x_max, P_apple_asY[Q_fit_asX == x_max], 
    labels = "D_apple", pos = 4)
#### save the graph into a file
# dev.copy(jpeg, file = "./outputs/graph_demand_curve_apple.jpg")
# dev.off()
### end #############################################
#### R script ########################################
# Last modfied: 2024-08-02
# File name: test_P_apple.R
# Keywords: 
#  test_P_apple: test the elastcity of P_apple,
#    using the demand function of Apple,
#
# Written by: Hiroshi Murao
#
# Purpose: Perform a hypothesis test  
#     for a model of the apple demand curve. 
#
# Original data source: ƌvNñ published by ӝvÁE
# The data can be obtained by visiting the following Web sites.
# ӝvǁAĉƌvN??Eҍ،ˉʁAēѐlȏÂ̐сAĕݐ̂PONxN??E݌v\ADƂPV\@
# H׿u̍؁ECєiжN̍؁jv`uʕȁv
#
# The following is a list of variable names in the data set. 
#
# year
# EXP_fresh_fruitQ_fresh_fruitPN_fresh_fruitP_fresh_fruit
# EXP_appleQ_applePN_appleP_apple
# EXP_mandarinQ_mandarinPN_mandarinP_mandarin
# EXP_pearQ_pearPN_pearP_pear
# EXP_grapeQ_grapePN_grapeP_grape
# EXP_persimmonQ_pwersimmonPN_persimmonP_persimmon
# EXP_peachQ_peachPN_peachP_peach
# EXP_watermelonQ_watermelonPN_watermelonP_watermelon
# EXP_melonQ_melonPN_melonP_melon
# EXP_strawberryQ_strawberryPN_strawberryP_strawberry
# EXP_bananaQ_bananaPN_bananaP_banana
# EXP_other_fruitQ_other_fruitPN_other_fruitP_other_fruit
# CPI
# income
# trend
#
#### first things to do ###############################
setwd("C:/book_SIMU/data_works/Apple/R_prg")
#### data preparations
in_data <- read.csv("./data/data_apple_JP.csv")
#print(cbind(in_data$year, in_data$Q_apple))
names(in_data)
print(in_data$Q_apple)
logQ_apple <- log(in_data$Q_apple)
logP_apple <- log(in_data$P_apple)
logIncome <- log(in_data$income)
logP_persimmon <- log(in_data$P_persimmon)
logP_mandarin <- log(in_data$P_mandarin)
trend <- in_data$trend
#### regression for H1, using lm()
model_H1 <- lm(logQ_apple ~ logP_apple + logIncome + logP_persimmon + trend)
summary(model_H1)
summary_H1 <- summary(model_H1)
names(summary_H1)
print(summary_H1$coefficients)
names(summary_H1$coefficients)
names(model_H1)
names(summary_H1)
num_df <- summary_H1$df
u <- summary_H1$residuals
# "sigma" is "Residual standard error".
SE_u <- summary_H1$sigma
# make sure that sigma_u = "Residual standard error"
summary(model_H1)
print(SE_u)
# SSR_H1 <- sum((logQ_apple - model_H1$fitted.values)^2)
SSR_H1 <- sum(u^2)
print(SSR_H1)
num_T <- length(logQ_apple)
print(num_T)
num_para <- 5
print(num_T - num_para)
print(num_df)
num_df <- num_df[2]
print(num_df)
#### t-test ####################################
# get SE_logP_apple
vcov_out <- vcov(model_H1)
print(vcov_out)
SE_logP_apple <- sqrt(vcov_out[2, 2])
print(SE_logP_apple)
# reference information
print(summary_H1$coefficients)
# SE_logP_apple = 0.140932
# get coef_logP_apple
print(model_H1$coefficients)
coef_logP_apple <- model_H1$coefficients[2]
print(coef_logP_apple)
#### get t-stat 
P_apple_abs <- abs(coef_logP_apple)
print(P_apple_abs)
t_stat <- (abs(coef_logP_apple) - 1)/SE_logP_apple
names(t_stat) <- "t_stat"
print(t_stat)
#### reference information 
print(summary_H1$coefficients)
## t critical value for a left-tailed test
t_critical <- qt(p = 0.05, df = num_df, lower.tail = TRUE) 
print("critical value for 5% significance") 
print(t_critical) 
## Comments on the t-test
# (1) t_stat is out side of t_critical value => Reject H0 and accept H1.
#### prepare for the F-test ####################
logQxP_apple <- logQ_apple + logP_apple
model_H0 <- lm(logQxP_apple ~ logIncome + logP_persimmon + trend)
u_H0 <- model_H0$residual
SSR_H0 <- sum(u_H0^2)
print(SSR_H0)
print(SSR_H1)
### F-test ####################################
J <- 1
numerator <- (SSR_H0 - SSR_H1)/J
denominator_1 <- SSR_H1/num_df
denominator_2 <- SE_u^2
# make sure that "denominator_1 = denominator_2".
print(denominator_1)
print(denominator_2)
denominator <- denominator_2
print(numerator)
F_stat <- numerator / denominator
names(F_stat) <- "F_stat"
print(F_stat)
## F critical value for a left-tailed test
F_critical <- qf(p = 0.05, df1 = 1, df2 = num_df, lower.tail = FALSE) 
print("critical value for 5% significance") 
print(F_critical) 
## Comments on the F-test
# (1) F_stat is out side of F_critical value => Reject H0 and accept H1.
#### compare the t-test and the F-test ##############
sq_t_stat <- t_stat^2
names(sq_t_stat) <- "sq_t_stat"
print(t_stat)
print(sq_t_stat)
print(F_stat)
#### end #############################################
### R Script ########################################
# Last modfied: 2025-07-03
# Name of script: graph_demand_curve_apple.R
# Keywords: 
# graph_demand_curve: get the graph of demand curve,
# apple: the demand curve of apples,
#
# Written by: Hiroshi Murao
#
# Purpose: regression analysis about apple demand in Japan. 
#
# Original data source: ƌvNñ published by ӝvÁE
# The data can be obtained by visiting the following Web sites.
# ӝvǁAĉƌvN??Eҍ،ˉʁAēѐlȏÂ̐сAĕݐ̂PONxN??E݌v\ADƂPV\@
# H׿u̍؁ECєiжN̍؁jv`uʕȁv
#
# The following is a list of variable names in the data set. 
#
# year
# EXP_fresh_fruitQ_fresh_fruitPN_fresh_fruitP_fresh_fruit
# EXP_appleQ_applePN_appleP_apple
# EXP_mandarinQ_mandarinPN_mandarinP_mandarin
# EXP_pearQ_pearPN_pearP_pear
# EXP_grapeQ_grapePN_grapeP_grape
# EXP_persimmonQ_pwersimmonPN_persimmonP_persimmon
# EXP_peachQ_peachPN_peachP_peach
# EXP_watermelonQ_watermelonPN_watermelonP_watermelon
# EXP_melonQ_melonPN_melonP_melon
# EXP_strawberryQ_strawberryPN_strawberryP_strawberry
# EXP_bananaQ_bananaPN_bananaP_banana
# EXP_other_fruitQ_other_fruitPN_other_fruitP_other_fruit
# CPI
# income
# trend
#
# Input: data of the above variables.
# Output: estimation results and estimated curves
#
##########################################################
### first things to do 
rm(list = ls())
setwd("C:/book_SIMU/data_works/Apple/R_prg")
### read data
in_data <- read.csv("./data/data_apple_JP.csv")
names(in_data)
# print(in_data)
### estimate the apple demand curve
out_lm = lm(log(Q_apple) ~ log(P_apple) + log(income) + log(P_persimmon) 
          + trend, data = in_data)
trend2pw <- in_data$trend^2
# out_lm = lm(log(Q_apple) ~ log(P_apple) + log(income) + log(P_persimmon) 
#            + trend + trend2pw, data = in_data)   
summary(out_lm)
#### preparing for plotting curves
b <- coefficients(out_lm)
resid <- residuals(out_lm)
# print(b)
logQ_fit_pA <- b[1] + 
               b[2]*(log(in_data$P_apple)) + 
               b[3]*mean(log(in_data$income)) +  
               b[4]*mean(log(in_data$P_persimmon) + 
               b[5]*mean(in_data$trend))
Q_fit_pA <- exp(logQ_fit_pA)
logQ_dot_pA <- logQ_fit_pA + resid
Q_dot_pA <- exp(logQ_dot_pA)
#### plot: logQ as y-axis, logP as x-axis, with matplot()
logP_apple <- log(in_data$P_apple)
print(logP_apple)
x_min <- min(logP_apple) 
x_max <- max(logP_apple)
print(cbind(x_min, x_max))
# 3.548647 3.862633
logQ_mix <- data.frame(logQ_fit_pA = logQ_fit_pA, 
                       logQ_dot_pA = logQ_dot_pA)
y_min <- min(logQ_mix$logQ_dot_pA) 
y_max <- max(logQ_mix$logQ_dot_pA)
print(cbind(y_min, y_max))
# 9.930534 10.19273
#### inverse demand curve
plot.new()
matplot(logP_apple, logQ_mix, 
  type = c("l", "p"), 
  lwd = 2, pch = 20, col = 1, 
  xlim = c(3.5, 3.9), 
  ylim = c(9.9, 10.2), 
  xlab = "log of apple price [yen/100g]", 
  ylab = "log of apple quantity [g]", 
  main = "inverse demand curve in log scale")
#### additional infromation
# text(x_max, logQ_mix$logQ_fit_aA[logP_apple == x_max], 
#    labels = "invD_apple", pos = 4)
#### save as a file
# dev.copy(jpeg, file = "./outputs/inv_demand_curve_logScale.jpg")
# dev.off()
#### plotting Q as y-axis, P as x-axis, with matplot()
xy_mix <- data.frame(P_apple = in_data$P_apple, 
                     Q_fit_pA = Q_fit_pA, 
                     Q_dot_pA = Q_dot_pA) 
print(xy_mix)
#### sort by P_apple
sortBy_pA <- xy_mix[order(xy_mix$P_apple), ]
x_min <- min(sortBy_pA$P_apple)
x_max <- max(sortBy_pA$P_apple)
print(cbind(x_min, x_max))
# 34.76626 47.59049
y_min <- min(sortBy_pA$Q_dot_pA)
y_max <- max(sortBy_pA$Q_dot_pA)
print(cbind(y_min, y_max))
# 20548.31 26708.34
print(sortBy_pA)
P_apple_asX <- data.frame(P_apple = sortBy_pA$P_apple)
print(P_apple_asX)
Q_mix_asY <- data.frame(Q_fit_aA = sortBy_pA$Q_fit_pA, 
                        Q_dot_pA = sortBy_pA$Q_dot_pA)
print(Q_mix_asY)
plot.new()
matplot(P_apple_asX, 
        Q_mix_asY, 
  type = c("l", "p"), 
  lwd = 2, pch = 20, col = 1, 
  xlim = c(34.0, 48.0), 
  ylim = c(20000, 28000), 
  xlab = "apple price [yen/100g]", 
  ylab = "apple quantity [g]", 
  main = "inverse per-capita demand curve in raw scale")
#### additional infromation
# text(x_max, Q_mix_asY$Q_fit_aA[P_apple_asX == x_max], 
#    labels = "invD_apple", pos = 4)
#### save the graph as a file
# dev.copy(jpeg, file = "outputs/inv_demand_curve_rawScale_wDot.jpg")
# dev.off()
#### plotting P as y-axis, Q as x-axis with plot()
#### sort by Q_apple
sortBy_qA <- xy_mix[order(xy_mix$Q_fit_pA), ]
Q_fit_asX <- sortBy_qA$Q_fit_pA
P_apple_asY <- sortBy_qA$P_apple
x_min <- min(Q_fit_asX)
x_max <- max(Q_fit_asX)
print(cbind(x_min, x_max))
# 20830.26 25118.47
y_max <- max(P_apple_asY)
y_min <- min(P_apple_asY)
print(cbind(y_min, y_max))
# 34.76626 47.59049
#### demand curve in raw scale
plot.new()
plot(Q_fit_asX, P_apple_asY, 
  type = "l", lwd = 2, col = 1, 
  xlim = c(20000, 27000), 
  ylim = c(34.0, 48.0), 
  xlab = "apple quantity [g]", 
  ylab = "apple price [yen/100g]", 
  main = "demand curve in raw scale")
#### additional information
text(x_max, P_apple_asY[Q_fit_asX == x_max], 
    labels = "D_apple", pos = 4)
#### save the graph into a file
# dev.copy(jpeg, file = "./outputs/graph_demand_curve_apple.jpg")
# dev.off()
### end #############################################
q()
