#### estimate VAR(p) model #########################
num_lag <- 6
out_VAR <- vars::VAR(USA[, c("x", "pi", "i")], p = num_lag, type = "const")
roots(out_VAR, modulus = TRUE)
# max modulus = 0.967835386 => OK.
summary(out_VAR)
#### preparing for SVAR(p) ###########################
# num_y = 3 => n(n-1)/2 = 3 restrictions
#### 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-10-07
# File name: graph_3_rho.R
# Keywords: 
#    graph_3_rho: graph of 3 integrated process,
#      rho = (0.8, 1.0, 1.2),
#
# Written by: Hiroshi Murao
#
#### first things to do 
setwd("C:/book_SIMU/data_works/Others/R_prg")
search()
#### make time data
time <- seq(-1, 10, by=1)
size_T <- length(time)
#### make y
y1 <- matrix(0, nrow = size_T, ncol = 1)
y2 <- matrix(0, nrow = size_T, ncol = 1)
y3 <- matrix(0, nrow = size_T, ncol = 1)
rho1 <- 0.8
rho2 <- 1.0
rho3 <- 1.2
cumY <- matrix(0, nrow = size_T, ncol = 1)
#### for-loop to compute I(2) process
u <- 1.0
y1[2] <- u
y2[2] <- u
y3[2] <- u
for (t in (3:size_T)){
 y1[t] <- rho1*y1[t-1] 
 y2[t] <- rho2*y2[t-1]
 y3[t] <- rho3*y3[t-1]
}
#### end for-loop
#### preparing for graph
print(y1)
y1_max <- max(y1)
print(y1_max)
x_data <- time[1:size_T]
y1_out <- y1[1:size_T]
y2_out <- y2[1:size_T]
y3_out <- y3[1:size_T]
y_data <- cbind(y1_out, y2_out, y3_out)
print(x_data)
print(y_data)
x_max <- max(x_data)
y1_max <- max(y_data[, 1])
y2_max <- max(y_data[, 2])
y3_max <- max(y_data[, 3])
print(c(y1_max, y2_max, y3_max))
#### plot graph
plot.new()
matplot(x_data, y_data, 
type = "l", lwd = 3, col = 1, 
lty = c("solid", "longdash", "dotted"),
xlim = c(-1, (size_T + 2)), 
ylim = c(-1, 7), 
xlab = "time", 
ylab = "marginal responses", 
main = "Three Types of Marginal Responses")
#### additional information to the graph
abline(h  =  0, col  =  "gray")
text(x_max, y_data[(x_data == x_max), 1], labels = "rho = 0.8", pos = 4)
text(x_max, y_data[(x_data == x_max), 2], labels = "rho = 1.0", pos = 4)
text(x_max, y_data[(x_data == x_max), 3], labels = "rho = 1.2", pos = 4)
#### save the graph as a file
# dev.copy(jpeg, file  =  "./outputs/three_rho.jpg")
# dev.off()
###### end #############################################
#### R script ############################################
# Update: 2023-10-09
# File name: graph_3_integrated_procss.R 
# Keywords: 
#    graph: graph of impulse response curves,  
#    3_integrated_procss: 3 I(d) process, d = (0, 1, 2),
#
# Written by: Hiroshi Murao
#
#### first things to do 
setwd("C:/book_SIMU/data_works/Others/R_prg")
search()
#### make time data
time <- seq(-3, 10, by = 1)
num_T <- length(time)
#### make the matrices for y's 
y0 <- matrix(0, nrow = num_T, ncol = 1)
y1 <- matrix(0, nrow = num_T, ncol = 1)
y2 <- matrix(0, nrow = num_T, ncol = 1)
y3 <- matrix(0, nrow = num_T, ncol = 1)
#### for-loop to compute impuluse response values
shock <- 1
y0[4] <- shock
y1[4] <- shock
y2[4] <- shock
y3[4] <- shock
for (t in (5:num_T)){
 y0[t] <- 0
 y1[t] <- y1[t-1]
 y2[t] <- 2*y2[t-1] - y2[t-2]
 y3[t] <- 3*y3[t-1] - 3*y3[t-2] + y3[t-3] 
}
#### end for-loop
#### make outputs 
num_T <- 10
out_y0 <- y0[3:num_T]
out_y1 <- y1[3:num_T]
out_y2 <- y2[3:num_T]
out_y3 <- y3[3:num_T]
# there are 8 values
print(out_y3)
#### prepare for a graph
x_data <- seq(-1, 6, by = 1)
y_data <- cbind(out_y0, out_y1, out_y2, out_y3)
print(x_data)
print(y_data)
x_max <- max(x_data)
y_max <- max(y_data[, 4])
y_min <- min(y_data[, 1])
print(c(y_min, y_max))
#### plot, using matplot()
plot.new()
matplot(x_data, y_data, 
type = "l", lwd = 2, col = 1, 
lty = c("dashed", "solid", "dashed", "solid"),
xlim = c(-1, (x_max+1)), 
ylim = c(-5, 30), 
xlab = "time", 
ylab = "impulse responses", 
main = "Impulse response curves")
# additional information to the graph
abline(h  =  0, col  =  "gray")
text(x_max, y_data[(x_data == x_max), 1], labels = "I(0)", pos = 4)
text(x_max, y_data[(x_data == x_max), 2], labels = "I(1)", pos = 3)
text(x_max, y_data[(x_data == x_max), 3], labels = "I(2)", pos = 4)
text(x_max, y_data[(x_data == x_max), 4], labels = "I(3)", 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/I0_I3_process.jpg")
# dev.off()
#### end #########################################################
q()
