q()
q()
q()
q()
q()
q()
q()
q()
q()
q()
q()
q()
q()
q()
q()
#### plot y_graph #################################
plot.new()
matplot(x_data, y_graph, 
    type = "l", lwd = 2, col = 1,
    lty = c("solid", "dotted", "dashed",  "dotdash"), 
    xlim = c(x_min, x_max + 3), 
    ylim = c(-0.20, 0.08), 
    xlab = "time (quarter)", 
    ylab = "impluse response of y1 (no unit)", 
    main = "Impluse response, SVAR Model, (P, Y, R, M)",
    sub = "shock = w3 = R, response = y1 = P" 
     )
## additional information to the graph
abline(h = 0, col = "gray")
text(x_max, y_graph[x_max], labels = "y1_s3", pos = 4)
# pos = 1: below,  pos = 2: left,  pos = 3: up,  pos = 4: right.
#### get y_graph 
# response = y1 = P, impluse = y3 = R => use of "y1_s3".
y_graph <- out_y[, 1] 
head(y_graph)
y_min <- min(y_graph)
y_max <- max(y_graph)
print(c(y_min, y_max))
#### get time index
num_H <- length(y_graph)
x_data <- seq(1, num_H, by = 1)
x_max <- max(x_data)
x_min <- min(x_data)
print(c(x_min, x_max))
#### plot y_graph #################################
plot.new()
matplot(x_data, y_graph, 
    type = "l", lwd = 2, col = 1,
    lty = c("solid", "dotted", "dashed",  "dotdash"), 
    xlim = c(x_min, x_max + 3), 
    ylim = c(-0.20, 0.08), 
    xlab = "time (quarter)", 
    ylab = "impluse response of y1 (no unit)", 
    main = "Impluse response, SVAR Model, (P, Y, R, M)",
    sub = "shock = w3 = R, response = y1 = P" 
     )
## additional information to the graph
abline(h = 0, col = "gray")
text(x_max, y_graph[x_max], labels = "y1_s3", pos = 4)
# pos = 1: below,  pos = 2: left,  pos = 3: up,  pos = 4: right.
q()
q()
setRepositories()
utils:::menuInstallPkgs()
devtools::install_github("FinYang/tsdl")
devtools::install_github("FinYang/tsdl")
library(tsdl)
library(prophet)
library(kableExtra)
#### Monthly Sales of Tasty Cola (Bowerman & O'Connell, 1993)
cola_monthly_sales <- subset(tsdl,"Sales")[6][[1]]
#### convert ts to df
start_date <- as.Date("1990-01-01") 
end_date <- length(cola_monthly_sales)
dates <- seq.Date(from = start_date, by = "month", length.out = end_date)
df <- data.frame(ds = dates, y = as.vector(cola_monthly_sales))
#### fit the model
model <- prophet(df)
#### make predictions
df_1993 <- make_future_dataframe(model, periods=12, freq='month')
predict_1993 <- predict(model, df_1993)
#### graphic presentations
plot(model, predict_1993)
install.packages("devtools")
devtools::install_github("FinYang/tsdl")
q()
q()
q()
q()
q()
q()
q()
q()
q()
q()
q()
q()
#### first things to do ###############################################
rm(list = ls())
setwd("c:/data_homework/OLS_mult/house_Kyoto/R_prg")
# setwd("c:/book_EMET/data_works/Used_car/R_prg")
#### read data
in_data <- read.csv("./data/data_house.csv")
names(in_data)
head(in_data)
#### make new variables
logP <- log(in_data$price)
landsize <- in_data$landsize
housesize <- in_data$housesize
age <- in_data$age
train <- in_data$train
bus <- in_data$bus
walk <- in_data$walk
dum_strcture <- in_data$structure
dum_garage <- in_data$garage
q()
q()
q()
q()
q()
q()
q()
q()
#### R script ########################################
# Last modfied: 2024-04-28
# File name: graph_line_dot.R
# Keywords for naming: 
#  graph_line_dot: make a graph with the estimated line and scatter dots, 
#     using a cobb-douglas production function,
#
# Written by: Hiroshi Murao
#
# Purpose: OLS estimation and related tasks.
# using the data of primary metal production, SIC 33.
#
# Original data source: Table 7.1 in Greene (2000)
# Greene, William (2000), Econometric Analysis, 4th edition
#
# Input: data located in subdirectory "data".
# Output: graphs located in subdirectory "outputs".
#
#######################################################
#### first things to do 
setwd("C:/book_SIMU/data_works/Metal/R_prg")
### read data
in_data <- read.csv("./data/data_metal_US.csv")
names(in_data)
head(in_data)
### OLS estimation
out_lm <- lm(log(output) ~ log(labor) + log(capital), data = in_data)
summary(out_lm)
names(out_lm)
summary_lm <- summary(out_lm)
names(summary_lm)
#### cheack some values
SSR_byF <- deviance(out_lm)
# deviance() function: the sum of deviance residuals squares is returned.
print(SSR_byF)
u <- summary_lm$residuals
resid <- residuals(out_lm)
num_T <- length(resid)
num_para <- 3
print(c(num_T, num_para))
sigma_byF <- SSR_byF/(num_T - num_para)
sigma_byM <- (t(resid)%*%resid)/(num_T - num_para)
print(sigma_byF)
print(sigma_byM)
# standard error of u
se_byF <- sqrt(sigma_byF)
print(se_byF)
SST <- deviance(lm(log(output) ~ 1, data = in_data))
print(SST)
Rsq <- 1 - SSR_byF/SST
print(Rsq)
summary(out_lm)
### preparations for plotting curves
b <- coefficients(out_lm)
print(b)
logY_labor <- b[1] + 
              b[2]*(log(in_data$labor)) + 
              b[3]*mean(log(in_data$capital))    
# fix capital at its mean
logY_labor_u <- logY_labor + resid
y_labor <- exp(logY_labor)
y_labor_u <- exp(logY_labor_u)
logY_capital <- b[1] + 
                b[2]*mean(log(in_data$labor)) + 
                b[3]*log(in_data$capital)     
# fix labor at its mean
logY_capital_u <- logY_capital + resid
y_capital <- exp(logY_capital)
y_capital_u <- exp(logY_capital_u)
#### plot (1): output vs. labor, with sort #################################
xy_labor <- data.frame(obs = in_data$obs, 
                       labor = in_data$labor, 
                       y_labor = y_labor, 
                       y_labor_u = y_labor_u)
print(xy_labor)
sortBy_labor <- xy_labor[order(xy_labor$labor), ]
print(sortBy_labor)
labor_min <- min(sortBy_labor$labor)
labor_max <- max(sortBy_labor$labor)
output_min <- min(sortBy_labor$y_labor_u)
output_max <- max(sortBy_labor$y_labor_u)
plot.new()
# with(sortBy_labor, 
# plot(labor, y_labor, 
# type = "l", lwd = 2, col = "red", 
# xlim = c(labor_min, labor_max), 
# ylim = c(output_min, output_max), 
# axes = F, xlab = "", ylab = ""))
#
# par(new = T)  # overwrite graphs
# with(sortBy_labor, 
# plot(labor, y_lab, 
# type = "p", pch = 20, col = "black", 
# xlim = c(labor_min, labor_max), 
# ylim = c(output_min, output_max), 
# xlab = "labor", 
# ylab = "output", 
# main = "output vs. labor where capital is controlled")
# ) # end of with()
#### the following line is equivatent to the above 3 commands.
# using matplot() function
with(sortBy_labor, 
  matplot(labor, cbind(y_labor, y_labor_u), 
  type = c("l", "p"), 
  lwd = 2, pch = 20, 
  col = c("black", "black"), 
  xlim = c(labor_min, labor_max), 
  ylim = c(output_min, output_max), 
  xlab = "labor", 
  ylab = "output", 
  main = "output vs. labor where capital is controlled")
) # end of with() 
#### save the graph as a file
# dev.copy(jpeg, file = "./outputs/output_vs_labor_with_resid.jpg")
# dev.off()
#### plot (2): output vs. capital, with sort #######################################
xy_capital <- data.frame(obs = in_data$obs, 
                         capital = in_data$capital, 
                         y_capital = y_capital, 
                         y_capital_u = y_capital_u)
head(xy_capital)
sortBy_capital <- xy_capital[order(xy_capital$capital), ]
head(sortBy_capital)
capital_min <- min(sortBy_capital$capital)
capital_max <- max(sortBy_capital$capital)
output_min <- min(sortBy_capital$y_capital_u)
output_max <- max(sortBy_capital$y_capital_u)
plot.new()
# with(sortBy_capital, 
# plot(capital, y_fit_capital, 
# type = "l", lwd = 2, col = "red", 
# xlim = c(capital_min, capital_max), 
# ylim = c(output_min, output_max), 
# axes = F, xlab = "", ylab = "")) 
#
# par(new=T) # overwrite graphs
# with(sortBy_capital, 
# plot(capital, y_dot_capital, 
# type = "p", pch = 20, col = "black", 
# xlim = c(capital_min, capital_max), 
# ylim = c(output_min, output_max), 
# xlab = "capital", 
# ylab = "output", 
# main = "output vs. capital where labor is controlled")
# ) # end of with()
#### the following line is equivatent to the above 3 commands.
# using matplot() function
with(sortBy_capital, 
  matplot(capital, cbind(y_capital, y_capital_u), 
  type = c("l", "p"), 
  lwd = 2, pch = 20, 
  col = c("black", "black"), 
  xlim = c(capital_min, capital_max), 
  ylim = c(output_min, output_max), 
  xlab = "capital", 
  ylab = "output", 
  main = "output vs. capital where labor is controlled") # end of matplot()
) # end of with()
#### save the graph as a file
# dev.copy(jpeg, file="./outputs/output_vs_capital_with_resid.jpg")
# dev.off()
#### plot (3): a wrong graph ##################################################
# make a wrong graph: fitted line "y_labor" with scatter plot of actual observations
labor_max <- max(in_data$labor)
labor_min <- min(in_data$labor)
output_max <- max(in_data$output)
output_min <- min(in_data$output)
plot.new()
with(sortBy_labor, 
  plot(labor, y_labor, 
  type = "l", lwd = 2, col = "red", 
  xlim = c(labor_min, labor_max), 
  ylim = c(output_min, output_max), 
  axes = F, xlab = "", ylab = "")
) # end of with()
par(new = T)  # overwrite graphs
with(in_data, 
  plot(labor, output, 
  type = "p", pch = 20, col = "black", 
  xlim = c(labor_min, labor_max), 
  ylim = c(output_min, output_max), 
  xlab = "labor", 
  ylab = "output", 
  main="wrong graph")
) # end of with()
#### save the graph as a file
# dev.copy(jpeg, file="./outputs/wrong_graph.jpg")
# dev.off()
### end ########################################
#### R Script ########################################
# Last modfied: 2025-04-02
# File name: test_F_constant_return.R
# Keywords for naming: 
#   test_F: F test, 
#   constant_return: test of constant return,
#      using a cobb-douglas production function,
#
# Written by: Hiroshi Murao
#
# Purpose: perform the F test.
# 
# Original data source: Table 7.1 in Greene (2000)
# Content od data: about American primary metals SIC 33.
# Greene, William (2000), Econometric Analysis, 4th edition
#
# colnames: "obs", "output", "labor", "capital".
#
# Input: data located in subdirectory "data".
# Output: information about F_test 
#
############################################################
#### first things to do 
rm(list = ls()) 
setwd("C:/book_SIMU/data_works/Metal/R_prg")
#### get data
in_data <- read.csv("./data/data_metal_US.csv")
head(in_data)
#### estimate the unrestriced model ##########################
out_lm <- lm(log(in_data$output) ~ log(in_data$labor) + log(in_data$capital), 
         data = in_data)
summary(out_lm)
#### get SSR_H1
names(out_lm)
num_df_resid <- out_lm$df.residual
summary_lm <- summary(out_lm)
names(summary_lm)
summary_lm$sigma
# sigma is SE of u, i.e., sigma^2 is MSE(u).
SE_u <- (summary(out_lm))$sigma
SSR_lm <- (SE_u^2)*num_df_resid
SSR_dv <- deviance(out_lm)
# deviance() function: the sum of deviance residuals squares is returned.
print(c(SSR_lm, SSR_dv))
# choose SSR_H1
SSR_H1 <- SSR_dv
print(SSR_H1)
num_obs <- nrow(in_data)
print(num_obs)
names(out_lm)
num_df_u <- out_lm$df.residual
print(num_df_u)
num_df_H1 <- num_df_u 
sigma_H1 <- SSR_H1/num_df_H1
print(sigma_H1)
#### estimate the restriced model ###################################
# lm_H0 <- lm((log(in_data$output) - log(in_data$labor)) ~ (log(in_data$capital) - log(in_data$labor)), data=in_data)
# The above lm_H0 is different from the below lm_H0. 
logYrL <- log(in_data$output) - log(in_data$labor)
logKrL <- log(in_data$capital) - log(in_data$labor)
# data_H0 <- data.frame(logYrL=logYrL, logKrL=logKrL)
lm_H0 <- lm(logYrL ~ logKrL)
summary(lm_H0)
#### get SSR_H0
SSR_H0 <- deviance(lm_H0)
print(SSR_H0)
num_J <- 1
#### perform the F test using scalars ##############################
F_stat <- ((SSR_H0 - SSR_H1)/num_J)/sigma_H1 
print("F statistic")
print(F_stat)
print(SSR_H0)
print(SSR_H1)
print(sigma_H1)
p_value <- pf(F_stat, num_J, num_df_H1, lower.tail = FALSE)
print("p value")
print(p_value)
#### estimation using matrices ########################################
y <- log(in_data[,2])
x1 <- rep(1, num_obs)
x2 <- log(in_data[,3])
x3 <- log(in_data[,4])
X <- cbind(x1, x2, x3)
print(X)
b <- solve(t(X)%*%X)%*%t(X)%*%y
print(b)
#### F test using matrices ############################################
R <- t(c(0,1,1))
print(R)
q <- 1
# XX <- t(X)%*%X 
XX_inv <- solve(t(X)%*%X)
# RXXR <- R%*%inv_XX%*%t(R)
RXXR_inv <- solve(R%*%XX_inv%*%t(R))
print(RXXR_inv)
SSR_diff <- t(R%*%b - q)%*%RXXR_inv%*%(R%*%b - q)
F_stat_mat <- (SSR_diff/num_J)/sigma_H1
print("F statistic")
print(F_stat_mat)
print(F_stat)
p_value_mat <- pf(F_stat_mat, num_J, num_df_H1, lower.tail = FALSE)
print("p value")
print(p_value_mat)
print(p_value)
#### Likelihood Ratio Test, start from reading data ########################
#### get data
# in_data <- read.csv("./data/metal_US.csv")
#### estimating the unrestricted model
model_H1_fit <- lm(log(in_data$output) ~ log(in_data$labor) + log(in_data$capital), 
                data = in_data)
summary(model_H1_fit)
#### get SSR_H1
SSR_H1 <- deviance(model_H1_fit)
num_obs <- nrow(in_data)
sigma_H1 <- SSR_H1/num_obs
#### estimating the restricted model
logYrL <- log(in_data$output) - log(in_data$labor)
logKrL <- log(in_data$capital) - log(in_data$labor)
model_H0_fit <- lm(logYrL ~ logKrL)
summary(model_H0_fit)
#### get SSR_H0
SSR_H0 <- deviance(model_H0_fit)
sigma_H0 <- SSR_H0/num_obs
#### compute the likelihood ratio (LR) test statistic #########
LR_stat <- num_obs*(log(sigma_H0) - log(sigma_H1))
print(LR_stat)
#### compute the p-value under the chi squared distribution
num_J <- 1
p_value_LR <- pchisq(LR_stat, num_J, lower.tail = FALSE)
print("p value, using LR test statistic")
print(p_value_LR)
#### end #############################################
#### R Script ########################################
# Last modfied: 2025-07-02
# Name of script: test_F_constant_return_byMat.R
# Keywords: 
# test_F: perform the F test,  
# constant_return: test the constant return, 
#    using Cobb-Douglas production function
# byMat: using matricies,
#
# Written by: Hiroshi Murao
#
# Purpose: This program performs an F test.
# 
# Original data source: Table 7.1 in Greene (2000)
# Content od data: about American primary metals SIC 33.
# Greene, William (2000), Econometric Analysis, 4th edition
#
# colnames are "obs", "output", "labor", "capital".
#
# Input: data located in subdirectory "data".
# Output: information about F_test 
#
#####################################################
#### first things to do 
rm(list = ls()) 
setwd("C:/book_SIMU/data_works/Metal/R_prg")
#### get data 
in_data <- read.csv("./data/data_metal_US.csv")
print(in_data)
#### estimation using lm() function ################
fit_H1 <- lm(log(in_data$output) ~ log(in_data$labor) + log(in_data$capital), 
            data=in_data)
summary(fit_H1)
#### estimation using matrices ######################
num_obs <- nrow(in_data)
print(num_obs)
y <- log(in_data[,2])
x1 <- rep(1, num_obs)
x2 <- log(in_data[,3])
x3 <- log(in_data[,4])
X <- cbind(x1, x2, x3)
b <- solve(t(X)%*%X)%*%t(X)%*%y
print(b)
### important numbers for the F test
num_para <- 3
num_df_H1 <- num_obs - num_para
num_J <- 1
#### F test using matrices #########################
resid <- y - X%*%b
SSR_H1 <- t(resid)%*%resid
s2_H1 <- SSR_H1/num_df_H1
R <- t(c(0,1,1))
print(R)
q <- 1
# XX <- t(X)%*%X 
XX_inv <- solve(t(X)%*%X)
# RXXR <- R%*%XX_inv%*%t(R)
RXXR_inv <- solve(R%*%XX_inv%*%t(R))
print(RXXR_inv)
SSR_diff <- t(R%*%b - q)%*%RXXR_inv%*%(R%*%b - q)
F_stat <- (SSR_diff/num_J)/s2_H1
print("F statistic")
print(F_stat)
p_value <- pf(F_stat, num_J, num_df_H1, lower.tail = FALSE)
print("p value")
print(p_value)
#### end #############################################
q()
