#### make Amat
Amat <- diag(3)
# Amat[1, 2:4] <- NA
Amat[2, 1] <- NA
Amat[3, 1:2] <- NA
print(Amat) 
#### make Bmat
Bmat <- diag(3)
print(Bmat)
#### estimate SVAR(p) ################################
# using SVAR() in "vars"
out_SVAR <- SVAR(out_VAR, estmethod = "direct",  
                 Amat = Amat, Bmat = Bmat,
                 max.iter = 100) 
summary(out_SVAR)
#### perform impulse response analysis ##############
# impluse = "i" = R_rate
num_H <- 24
irf_SVAR <- irf(out_SVAR, 
                     impluse = "i", 
                     n.ahead = num_H, 
                     ortho = FALSE)
plot(irf_SVAR)
names(irf_SVAR)
names(irf_SVAR$irf)
head(irf_SVAR$irf$"i")
irf_s3 <- irf_SVAR$irf$"i"
head(irf_s3)
#### get irf-values, given "i" shock 
y1_s3 <- irf_s3[, 1]
y2_s3 <- irf_s3[, 2]
y3_s3 <- irf_s3[, 3]
head(y3_s3)
#### prepare for plotting ############
# get ys_s3
ys_s3 <- cbind(y1_s3, y2_s3, y3_s3)
# ys_s3 <- irf_s3
head(ys_s3)
ini_0 <- c(0, 0, 0)
ys_w_ini <- rbind(ini_0, ys_s3)
head(ys_w_ini)
colnames(ys_w_ini) <- c("y1_s3", "y2_s3", "y3_s3")
y_min <- min(ys_w_ini)
y_max <- max(ys_w_ini)
print(c(y_min, y_max))
print(nrow(ys_w_ini))
# get x_axis
num_row <- nrow(ys_w_ini)
x_axis <- matrix(-1:(num_row-2), nrow = num_row, ncol = 1)
x_min <- min(x_axis)
x_max <- max(x_axis)
print(c(x_min, x_max))
print(length(x_axis))
#### plot an x-y graph 
plot.new()
matplot(x_axis, ys_w_ini, 
type = "l", lwd = 2, col = 1, 
lty = c("solid", "dotdash", "dotted"), 
xlim = c(-1, (x_max+6)), 
ylim = c(y_min, y_max), 
xlab = "time (quarter)", 
ylab = "impulse responses (percent)", 
main = "SVAR, a positive R shock")
#### additional information to the graph
abline(h = 0, col = "gray")
text(x_max, ys_w_ini[x_axis == x_max, 1], labels = "Y_gap", pos = 4)
text(x_max, ys_w_ini[x_axis == x_max, 2], labels = "P_growth", pos = 4)
text(x_max, ys_w_ini[x_axis == x_max, 3], labels = "R", pos = 4)
# pos=1: below,  pos = 2: left,  pos = 3: up,  pos = 4: right.
#### save the graph as a file
# dev.copy(jpeg, file = "./outputs/irf_vars_SVAR_shock_R.jpg")
dev.copy(jpeg, file = "C:/book_SVAR/data_works/Harwartz/outputs/irf_vars_SVAR_shock_R.jpg")
dev.off()
# dev.copy(pdf, file = "./outputs/irf_var_SVAR_shock_R.pdf")
# dev.copy(postscript, file = "./outputs/irf_var_SVAR_shock_R.eps")
#### end #############################################
#### R script ########################################
# Last modfied: 2024-10-19
# File name: irf_var_SVEC_shock_R.R
# Keywords: 
#   irf: impulse response analysis, 
#   var_SVEC:  using vars::SVEC() function,  
#   shock_R: shock to R_rate, 
#      tighting monetary policy, 
#      using y = c("x", "pi", "i") = (Y_gap, P_growth, R_rate)  
#            
# Written by: Hiroshi Murao
#
# Purpose: perfome a comparion of impulse response analysis,
#    based on 3 models: a level-VAR model, a difference-VAR, and a VEC model. 
#    We transform difference-VAR(p-1) to level-VAR(p) after model estimation.
#    We transform VEC(p-1) to level-VAR(p) after the model estimatiion.  
#
# Data source:
#    USA in Package "svars",  
#
# Data information: 
#   Data type: U.S. quarterly data 
#   Time periods of original data set: 1965Q1 - 2008Q3
# 
# List of variables:
# x:  Log-deviation of real gross domestic product (GDP) from the potential output.     
# pi: Quarter-on-quarter growth rate of GDP deflator.    
# i: Federal funds rate.
#
#############################################
#### first things to do 
rm(list = ls())
setwd("C:/book_SVAR/data_works/Harwartz/R_prg")
library(urca) # unit root, cointegration
library(vars) # VAR
library(svars) # SVAR 
search()
#### get data
data("USA")
head(USA)
plot(USA) 
#### make new variables
y1 <- USA[, "x"]
y2 <- USA[, "pi"]
y3 <- USA[, "i"]
head(cbind(y1, y2, y3))
#### select num_lag
out_VARselect <- vars::VARselect(USA[, c("x", "pi", "i")], lag.max = 10, type = "const")
print(out_VARselect)
# AIC(n)  HQ(n)  SC(n) FPE(n) 
#     6      3      3      6 
#### estimate SVEC model ################################
## prepare for the estimation of VEC(p), using ur.jo() in "urca"
num_lag <- 6
out_cajo <- ca.jo(x = USA[, c("x", "pi", "i")], type = "eigen", 
                  ecdet = "const", K = num_lag, spec = "transitory")
summary(out_cajo)
## make LRmat
LRmat <- matrix(NA, nrow = 3, ncol = 3)
LRmat[1:3, 3] <- 0
LRmat[1, 2:3] <- 0
print(LRmat) 
## make SLmat
SRmat <- matrix(NA, nrow = 3, ncol = 3)
print(SRmat)
out_SVEC <- vars::SVEC(out_cajo, r = 1, LR = LRmat, SR = SRmat, 
                  lrtest = FALSE)
summary(out_SVEC)
#### impluse responce analysis, using irf() function in Package "vars"
num_H <- 24 
irf_SVEC <- irf(out_SVEC, impulse = "i", n.ahead = num_H, 
               ortho = FALSE, cumulative = FALSE, 
               boot = TRUE, ci = 0.95, runs = 500)
plot(irf_SVEC) 
#### get a better graph 
# checking how to pick up some outputs
# str(summary(irf_SVEC)) 
# irf of all_Y from "i"
print(irf_SVEC$"irf"$"i") 
irf_s3 <- irf_SVEC$"irf"$"i"
head(irf_s3)
#### get irf-values, given "i" shock 
y1_s3 <- irf_s3[, 1]
y2_s3 <- irf_s3[, 2]
y3_s3 <- irf_s3[, 3]
head(y3_s3)
#### prepare for plotting ############
# get ys_s3
ys_s3 <- cbind(y1_s3, y2_s3, y3_s3)
# ys_s3 <- irf_s3
head(ys_s3)
ini_0 <- c(0, 0, 0)
ys_w_ini <- rbind(ini_0, ys_s3)
head(ys_w_ini)
colnames(ys_w_ini) <- c("y1_s3", "y2_s3", "y3_s3")
y_min <- min(ys_w_ini)
y_max <- max(ys_w_ini)
print(c(y_min, y_max))
print(nrow(ys_w_ini))
# get x_axis
num_row <- nrow(ys_w_ini)
x_axis <- matrix(-1:(num_row-2), nrow = num_row, ncol = 1)
x_min <- min(x_axis)
x_max <- max(x_axis)
print(c(x_min, x_max))
print(length(x_axis))
#### plot an x-y graph 
plot.new()
matplot(x_axis, ys_w_ini, 
type = "l", lwd = 2, col = 1, 
lty = c("solid", "dotdash", "dotted"), 
xlim = c(-1, (x_max+6)), 
ylim = c(y_min, y_max), 
xlab = "time (quarter)", 
ylab = "impulse responses (percent)", 
main = "SVEC, a positive R shock")
#### additional information to the graph
abline(h = 0, col = "gray")
text(x_max, ys_w_ini[x_axis == x_max, 1], labels = "Y_gap", pos = 4)
text(x_max, ys_w_ini[x_axis == x_max, 2], labels = "P_growth", pos = 4)
text(x_max, ys_w_ini[x_axis == x_max, 3], labels = "R", pos = 3)
# pos=1: below,  pos = 2: left,  pos = 3: up,  pos = 4: right.
#### save the graph as a file
# dev.copy(jpeg, file = "./outputs/irf_var_SVEC_shock_R.jpg")
dev.copy(jpeg, file = "C:/book_SVAR/data_works/Harwartz/outputs/irf_vars_SVEC_shock_R.jpg")
dev.off()
# dev.copy(pdf, file = "./outputs/irf_var_SVEC_shock_R.pdf")
# dev.copy(postscript, file = "./outputs/irf_var_SVEC_shock_R.eps") 
#### end #############################################
q()
q()
q()
q()
q()
q()
q()
q()
q()
q()
q()
q()
q()
q()
q()
q()
q()
q()
q()
q()
q()
q()
q()
q()
q()
q()
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()
q()
q()
q()
q()
q()
q()
q()
q()
q()
#### R script ########################################
# Last modfied: 2023-09-13
# File name: grph_line_dot_wError.R
# Keywords: 
#   graph: get an illustrative figure for regression,
#   line_dot: line with dots,
#   wError: error will be addted to the graoh by hands,     
#            
# Written by: Hiroshi Murao
#
# Purpose: perfome the simple regression,
#      using one of equations useded in the Klein Model.  
#
# Data source:
#  Klein's (1950) Model I of the US economy.
#  Klein, L. (1950). Economic Fluctuations in the United States, 
#  1921-1941. John Wiley, New York.
#
# List of variables: names used in bimets Package.
# (1) cn as Private Consumption Expenditure;
# (2) i as Investment;
# (3) w1 as Wage Bill of the Private Sector (Demand for Labor);
# (4) p as Profits;
# (5) k as Stock of Capital Goods;
# (6) y as Gross National Product;
# (7) w2 as Wage Bill of the Government Sector;
# (8) time as an annual index of the passage of time;
# (9) g as Government Expenditure plus Net Exports;
# (10) t as Business Taxes.
#
##########################################################
#### first things to do 
rm(list = ls())
setwd("C:/book_SIMU/data_works/Klein/R_prg")
library(bimets) # Time Series and Econometric Modeling
library(systemfit) # Time Series and Econometric Modeling
search()
#### get the data for Klein's (1950) Model I of the US economy.
in_data <- read.csv("./data/data_Klein_bimets.csv", 
           na.strings = c("", ".", "NA", "#N/A"), 
           stringsAsFactors = FALSE,
           colClasses = c(rep("numeric", 10)))
head(in_data)
#### make variables for use
y <- in_data$cn
x <- in_data$y
#### scatter plots plus the regression line
y_min <- min(y)
y_max <- max(y)
print(c(y_min, y_max))
# [1] 39.8 69.7
x_min <- min(x)
x_max <- max(x)
print(c(x_min, x_max))
# [1] 40.6 85.3
# scatter plot 
plot.new()
matplot(x, y, 
type = "p", pch = 16, col = 1,  
xlim = c(30.0, 90.0), 
ylim = c(30.0, 80.00), 
xlab = "x", 
ylab = "y")
### additional information to the graph
abline(lm(y ~ x))
# text(x_max, u3_asY[(u1_hat_asX == x_max), 1], labels = "Reg_line", pos = 3)
# pos = 1: below,  pos = 2: left,  pos = 3: up,  pos = 4: right.
#### save the graph into a file
# dev.copy(jpeg, file = "./outputs/reg_illustrative_fig_pch16.jpg")
# dev.off()
#### Comments
# (1) 
# (2)
#### end #############################################
#### R script ########################################
# Last modfied: 2023-09-13
# File name: grph_line_dot_wError.R
# Keywords: 
#   graph: get an illustrative figure showing a regression line,
#   line_dot: line with dots, dats = data points,
#   wError: an error line (a vertical line) will be added to the graph by hands,     
#            
# Written by: Hiroshi Murao
#
# Purpose: perfome the simple regression,
#      using one of equations useded in the Klein Model.  
#
# Data source:
#  Klein's (1950) Model I of the US economy.
#  Klein, L. (1950). Economic Fluctuations in the United States, 
#  1921-1941. John Wiley, New York.
#
# List of variables: names used in bimets Package.
# (1) cn as Private Consumption Expenditure;
# (2) i as Investment;
# (3) w1 as Wage Bill of the Private Sector (Demand for Labor);
# (4) p as Profits;
# (5) k as Stock of Capital Goods;
# (6) y as Gross National Product;
# (7) w2 as Wage Bill of the Government Sector;
# (8) time as an annual index of the passage of time;
# (9) g as Government Expenditure plus Net Exports;
# (10) t as Business Taxes.
#
##########################################################
#### first things to do 
rm(list = ls())
setwd("C:/book_SIMU/data_works/Klein/R_prg")
library(bimets) # Time Series and Econometric Modeling
library(systemfit) # Time Series and Econometric Modeling
search()
#### get the data for Klein's (1950) Model I of the US economy.
in_data <- read.csv("./data/data_Klein_bimets.csv", 
           na.strings = c("", ".", "NA", "#N/A"), 
           stringsAsFactors = FALSE,
           colClasses = c(rep("numeric", 10)))
head(in_data)
#### make variables for use
y <- in_data$cn
x <- in_data$y
#### scatter plots plus the regression line
y_min <- min(y)
y_max <- max(y)
print(c(y_min, y_max))
# [1] 39.8 69.7
x_min <- min(x)
x_max <- max(x)
print(c(x_min, x_max))
# [1] 40.6 85.3
# scatter plot 
plot.new()
matplot(x, y, 
type = "p", pch = 16, col = 1,  
xlim = c(30.0, 90.0), 
ylim = c(30.0, 80.00), 
xlab = "x", 
ylab = "y")
### additional information to the graph
abline(lm(y ~ x))
# text(x_max, u3_asY[(u1_hat_asX == x_max), 1], labels = "Reg_line", pos = 3)
# pos = 1: below,  pos = 2: left,  pos = 3: up,  pos = 4: right.
#### save the graph into a file
# dev.copy(jpeg, file = "./outputs/graph_line_dot_wError.jpg")
# dev.off()
#### Comments
# (1) 
# (2)
#### end #############################################
q()
