Y_diff <- na.omit(diff(Y, lag = 1))
head(Y)
head(Y_lag)
CP_lag <- CP_lag[-length(CP_lag)]
IP_lag <- IP_lag[-length(IP_lag)]
R_lag <- R_lag[-length(R_lag)]
MrP_lag <- MrP_lag[-length(MrP_lag)]
Y_lag <- Y_lag[-length(Y_lag)]
# We deal with data-length problem later on.  
#### get seasonal dummies for "bimets"
dum_Q1 <- as.numeric(data_use$dum_Q1)
dum_Q2 <- as.numeric(data_use$dum_Q2)
dum_Q3 <- as.numeric(data_use$dum_Q3)
#### check data sizes for #bimets"
print(c(length(date), length(year), length(quarter), 
      length(CP), length(IP), length(R), length(Y), 
      length(MrP), length(G), length(TAX), length(YD),
      length(dum_Q1), 
      length(Y_lag), length(Y_diff)))
#### check the contents of some variables 
plot(R, type = "l", lwd = 2)
plot(TAX, type = "l", lwd = 2)
plot(G, type = "l", lwd = 2)
#### set data_obs for "bimets" ###################################
# date <- date[-1] # 'closure' ^̃IuWFNg͕ԕʑÓ܉\ł܂قсB͂
data_obs <- data.frame(date = date[-1], year = year[-1], quarter = quarter[-1],   
            CP = CP[-1], IP = IP[-1], R = R[-1], Y = Y[-1],
            MrP = MrP[-1], G = G[-1], TAX = TAX[-1], YD = YD[-1],
            dum_Q1 = dum_Q1[-1], dum_Q2 = dum_Q2[-1], dum_Q3 = dum_Q3[-1],
            CP_lag = CP_lag, IP_lag = IP_lag, R_lag = R_lag, 
            MrP_lag = MrP_lag, Y_lag = Y_lag
            )
# we use CP_lag in "data_obs" for lm() function in different R scripts as well.
with(data_obs, 
print(c(length(date), length(year), length(quarter), 
      length(CP), length(IP), length(R), length(Y), 
      length(MrP), length(G), length(TAX), length(YD),
      length(dum_Q1), 
      length(CP_lag), length(IP_lag), length(MrP_lag)
)))
# check lagged values with no lagged values
head(data_obs$Y)
head(data_obs$Y_lag)
# check "data_obs"
head(data_obs) # starting 2003Q3 
tail(data_obs) # ending 2022Q3
#### OLS estimation ########################
OLS_out <- lm(R ~ Y + MrP + R_lag,
              data = data_obs)
# We don't use seasonal dummies for "R" equation.
summary(OLS_out) 
names(OLS_out)
# coefficient of "Y"
beta_OLS <- OLS_out$coefficients[2]
print(beta_OLS)
#### TSLS estimation ########################
Y_est <- lm(Y ~ G + TAX + MrP_lag
            + CP_lag + IP_lag + R_lag
            + dum_Q1 + dum_Q2 + dum_Q3, 
            data = data_obs)
summary(Y_est) 
Y_hat <- Y_est$fitted.values
head(Y_hat)
print(c(length(Y), length(Y_hat)))
data_TSLS <- data.frame(R = data_obs$R,
                        Y_hat = Y_hat,
                        MrP = data_obs$MrP, 
                        R_lag = data_obs$R_lag,
                        dum_Q1 = data_obs$dum_Q1, 
                        dum_Q2 = data_obs$dum_Q2, 
                        dum_Q3 = data_obs$dum_Q3
                        )
TSLS_out <- lm(R ~ Y_hat + MrP + R_lag,  
               data = data_TSLS)
# We don't use seasonal dummies for "R" equation.
summary(TSLS_out) 
names(TSLS_out) 
print(TSLS_out$coefficients)
# coefficient of "Y_hat"
beta_TSLS <- TSLS_out$coefficients[2]
print(beta_TSLS)
print(beta_OLS)
#### LIML estimation ########################
# Some inputs of ivmodel() function are shown:
# Y: A numeric vector of outcomes.
# D: A vector of endogenous variables.
# Z: A matrix or data frame of instruments.
# X: A matrix or data frame of (exogenous) covariates.
#
# The model is: Y = X*alpha + D*beta + error
# E(error given X and Z) = 0
# Main interest is: beta.
# (1) dependent variable of "R" equation
Y_ivm = data_obs[, "R"]
# names(Y_ivm) <- "R" 
head(Y_ivm)
# (2) one endogenous explanatory variable in "R" equation,
D_ivm <- data_obs[, "Y"]
# names(D_ivm) <- "Y" 
head(D_ivm)
# (3) intruments for D_ivm = Y
# one intrument case
# Z_ivm <- data_obs[, "Y_lag"]
# names(Z_ivm) <- "Y_lag" 
# more than two instruments case
# The choice of "Z_ivm" is important for the output.
# Z_ivm <- data_obs[, c("Y_lag", "G", "TAX")]
# colnames(Z_ivm) <- c("Y_lag", "G", "TAX")
# Z_ivm <- data_obs[, c("CP", "IP", "G")]
Z_ivm <- data_obs[, c("CP_lag", "IP_lag", "R_lag", "G", "TAX", "dum_Q1", "dum_Q2", "dum_Q3")]
head(Z_ivm)
# (4) exogenous explanatory variables in "R" equation,
X_name <- c("MrP", "R_lag")
X_ivm <- data_obs[, X_name]
# We don't use seasonal dummies for "R" equation.
# dum_Q1 = data_obs$dum_Q1, dum_Q2 = data_obs$dum_Q2, dum_Q3 = data_obs$dum_Q3
head(X_ivm)
#### input data to ivmodel() ###############################
head(cbind(Y_ivm, D_ivm, Z_ivm, X_ivm))
#### make "ivmodel" #########################################
out_ivm <- ivmodel(Y = Y_ivm, D = D_ivm, Z = Z_ivm, X = X_ivm, 
                   intercept = TRUE,
                   beta0 = 0, alpha = 0.05,
                   k = c(0, 1),
                   manyweakSE = TRUE, heteroSE = FALSE, clusterID = NULL,
                   deltarange = NULL, na.action = na.omit
                   ) 
#### important note #########################################
# (1) use "manyweakSE = TRUE"
# cheack the output of ivmodel()
print(out_ivm)
#### do LIML estimation #####################################
# The LIML estimator has a certain value of k.
out_LIML_ivm <- LIML(out_ivm,
                     beta0 = 0, alpha = 0.05,
                     manyweakSE = TRUE, heteroSE = FALSE,clusterID = NULL
                     )
# summary() does not work as expected.
summary(out_LIML_ivm) 
summary(out_LIML_ivm$point.est) 
print(out_LIML_ivm)
print(out_LIML_ivm$point.est)
print(out_LIML_ivm$point.est.other)
coef_out <- coef.ivmodel(out_LIML_ivm)
print(coef_out)
beta_LIML_ivm <- out_LIML_ivm$point.est
names(beta_LIML_ivm) <- "Y"
print(beta_LIML_ivm)
#### comparison of the coeffients of "Y" #####################
print(c(beta_OLS, beta_TSLS, beta_LIML_ivm[1]))
#### Comments ################################################
# (1) ivmodel() function is very sensitive to the choice of instruments "Z".
# (2) It may be better to use TSLS estiomation. 
#### end #####################################################
#### input data to ivmodel() ###############################
# Y_ivm: "MrP"
# D_ivm: "Y"
head(cbind(Y_ivm, D_ivm, Z_ivm, X_ivm))
#### R script ########################################
# Last modfied: 2025-08-16
# File name: est_ivmodel_LIML_Eq_M_noGood.R
# Keywords: 
#   est_ivmodel_LIML: estimate an equation, using LIML method,     
#        using ivmodel::LIML() function,
#   Eq_M: Money demand equation in an IS-LM model,
#   noGood: ivmodel::ivmodel() function might not work well,
#            
# Written by: Hiroshi Murao
#
# Purpose: perfomes the LIML estimation,
#      using ivmodel::LIML() function.  
#
# Data source:
# (1) Website of Japan Cabinet Office: 
#       https://www.esri.cao.go.jp/jp/sna/data/data_list/kakuhou/files/files_kakuhou.html
# (2) Website of the Bank of Japan: 
#       https://www.stat-search.boj.or.jp/ 
#
# Data information: 
#   Data type: Japan quarterly data 
#   Time periods of original data set: 1994:01 - 2022:03
#   Time periods of data set for use: 2003:02 - 2022:03
# 
# List of variables:
# (1) date : date information such as "1994Q1".
# (2) year : 
# (3) quarter : one of (1, 2, 3, 4)
# (4) GDP : Gross Domestic Product, unit: 10 billion yen 
# (5) C_prv : Consuption of private sector
# (6) I_prv_jutaku
# (7) I_prv_factory
# (8) I_prv = I_prv_jutaku + I_prv_factory : Investment of private sector
# (9) C_gov: Consuption of public sector
# (10) I_gov : Investment of public sector
# (11) Exports
# (12) Imports
# (13) GDI : Gross Domestic Income
# (14) GNI : Gross National Income
# (15) I_fixed = I_prv + I_gov = I_gross - I_invento
# (16) P_GDP : Deflator of GDP
# (17) P_C_prv : Deflator for C_prv
# (18) P_exports
# (19) P_imports
# (20) P_GNI : Deflator of GNI
# (21) M2 : Money stock M2, unit: 10 billion yen
# (22) M3
# (23) M_base: Monetary Base, unit: 10 billion yen
# (24) R_discount: Interest rate, discount rate, unit: percent (%),
# (25) R_call : Interest rate, call rate
# (26) R_loan_sogo : Interest rate, loan rate, "sogo" type
# (27) R_loan_short
# (28) R_loan_long
# (29) NDI : National Disporsable Income, unit: 10 billion yen
# (30) CPI: Consumer Price Index
# (31) IIP_sAjst: Index of Industrical Production, seasonal adjusted
#
##########################################################
#### first things to do 
rm(list = ls())
setwd("C:/book_SIMU/data_works/ISLM/R_prg")
library(ivmodel) # Fitting Intrument Variable Model
search()
#### get data
in_data <- read.csv("./data/data_ISLM.csv", 
           na.strings = c("", ".", "NA", "#N/A"), 
           stringsAsFactors = FALSE,
           colClasses = c("character", rep("numeric", 28)))
head(in_data)
sapply(in_data, class)
#### make seasonal dummies
num_inObs <- nrow(in_data)
print(num_inObs)
dum_Qs <- matrix(0, nrow = num_inObs, ncol = 4)
#### start t-loop
for (t in (1:num_inObs)){
  if (in_data[t, "quarter"] == 1) {
    dum_Qs[t, 1] <- 1
  } else if (in_data[t, "quarter"] == 2) {
    dum_Qs[t, 2] <- 1
  } else if (in_data[t, "quarter"] == 3) {
    dum_Qs[t, 3] <- 1
  } else dum_Qs[t, 4] <- 1
} 
#### end of t-loop
#### check dum_Qs
head(dum_Qs)
colnames(dum_Qs) <- c("dum_Q1", "dum_Q2", "dum_Q3", "dum_Q4")
dum_Qs <- dum_Qs[, -4]
colnames(dum_Qs)
head(dum_Qs)
#### marge dum_Qs to in_data
in_data <- cbind(as.matrix(in_data), as.matrix(dum_Qs))
in_data <- as.data.frame(in_data)
head(in_data) 
## choose the time periods without missing values
num_del <- 37
# num_del <- 57
data_noNA <- in_data[-(1:num_del), ]  
head(data_noNA)
# make sure that data_noNA starts from date = 2003Q2.
#### get data for use
data_use <- data_noNA
# data_use <- as.data.frame(subset(data_noNA, data_noNA$year >= 2004))
head(data_use)
#### get variables for possible use
# use as.nemeric() function if "colClasses" is not used.
# unit: trillion yen
# X*4 => growth rates would be annual rates
C_prv <- as.numeric(data_use$C_prv)/1000
I_prv <- as.numeric(data_use$I_prv)/1000
GDP <- as.numeric(data_use$GDP)/1000
IIP <- as.numeric(data_use$IIP_sAjst)/1000
M2 <- as.numeric(data_use$M2)/1000
M3 <- as.numeric(data_use$M3)/1000
M_base <- as.numeric(data_use$M_base)/1000
P_GDP <- as.numeric(data_use$P_GDP)
CPI <- as.numeric(data_use$CPI)
NDI <- as.numeric(data_use$NDI)/1000
# interest rates are already annual rates
R_discount <- as.numeric(data_use$R_discount) 
R_call <- as.numeric(data_use$R_call)
R_loan_sogo <- as.numeric(data_use$R_loan_sogo)
R_loan_short <- as.numeric(data_use$R_loan_short)
# check the content of interest rates
plot(R_call, type = "l", pch = 2, col = 1)
plot(R_loan_sogo, type = "l", pch = 2, col = 1)
plot(R_loan_short, type = "l", pch = 2, col = 1)
# important info.
# "R_call" has many zeroes with less variations.
# "R_loan_short" has more variations.
# seasonal dummy variables
dum_Qs <- data_use[, c("dum_Q1", "dum_Q2", "dum_Q3")]
head(dum_Qs)
# variables for time information
year_org <- as.numeric(data_use$year) 
quarter_org <- as.numeric(data_use$quarter)
#### convert date_info to date_asDate
date_info <- data_use$date
# date_info contains the value of form "1994Q2".
head(date_info)
## The following function is found from Internet.
date_asDate <- as.Date(
  sapply(strsplit(date_info, 'Q'), 
  function(x) paste(1, seq(1,10,3)[as.integer(x[2])], x[1], sep = '-')),
  format = '%d-%m-%Y')
## end of the function
#### check the output of as.Date() function
head(date_asDate)
# date with the format of "year-month-day", for example, "1994-04-01".
#### choose variables for the model #######################
year <- year_org
quarter <- quarter_org
date <- date_asDate
CP <- C_prv
IP <- I_prv
# R <- R_call
R <- R_loan_short
Y <- GDP
M <- M2
P <- P_GDP/100
MrP <- M/P
YD <- NDI
TAX <- Y - YD
G_id <- Y - CP - IP
CG <- as.numeric(data_use$C_gov)/1000
IG <- as.numeric(data_use$I_gov)/1000
EXPS <- as.numeric(data_use$Exports)/1000
IMPS <- as.numeric(data_use$Imports)/1000
G_noID <- CG + IG + (EXPS - IMPS)
# choose G between G_id and G_noID
# G <- G_id # We use the identity: Y = C + I + G
G <- G_noID
## lag and diff
CP_lag <- na.omit(lag(CP, lag = 1))
IP_lag <- na.omit(lag(IP, lag = 1))
R_lag <- na.omit(lag(R, lag = 1))
MrP_lag <- na.omit(lag(MrP, lag = 1))
Y_lag <- na.omit(lag(Y, lag = 1))
Y_diff <- na.omit(diff(Y, lag = 1))
head(Y)
head(Y_lag)
CP_lag <- CP_lag[-length(CP_lag)]
IP_lag <- IP_lag[-length(IP_lag)]
R_lag <- R_lag[-length(R_lag)]
MrP_lag <- MrP_lag[-length(MrP_lag)]
Y_lag <- Y_lag[-length(Y_lag)]
# We deal with data-length problem later on.  
#### get seasonal dummies for "bimets"
dum_Q1 <- as.numeric(data_use$dum_Q1)
dum_Q2 <- as.numeric(data_use$dum_Q2)
dum_Q3 <- as.numeric(data_use$dum_Q3)
#### check data sizes for #bimets"
print(c(length(date), length(year), length(quarter), 
      length(CP), length(IP), length(R), length(Y), 
      length(MrP), length(G), length(TAX), length(YD),
      length(dum_Q1), 
      length(Y_lag), length(Y_diff)))
#### check the contents of some variables 
plot(R, type = "l", lwd = 2)
plot(TAX, type = "l", lwd = 2)
plot(G, type = "l", lwd = 2)
#### set data_obs for "bimets" ###################################
# date <- date[-1] # 'closure' ^̃IuWFNg͕ԕʑÓ܉\ł܂قсB͂
data_obs <- data.frame(date = date[-1], year = year[-1], quarter = quarter[-1],   
            CP = CP[-1], IP = IP[-1], R = R[-1], Y = Y[-1],
            MrP = MrP[-1], G = G[-1], TAX = TAX[-1], YD = YD[-1],
            dum_Q1 = dum_Q1[-1], dum_Q2 = dum_Q2[-1], dum_Q3 = dum_Q3[-1],
            CP_lag = CP_lag, IP_lag = IP_lag, R_lag = R_lag, 
            MrP_lag = MrP_lag, Y_lag = Y_lag
            )
# we use CP_lag in "data_obs" for lm() function in different R scripts as well.
with(data_obs, 
print(c(length(date), length(year), length(quarter), 
      length(CP), length(IP), length(R), length(Y), 
      length(MrP), length(G), length(TAX), length(YD),
      length(dum_Q1), 
      length(CP_lag), length(IP_lag), length(MrP_lag)
)))
# check lagged values with no lagged values
head(data_obs$Y)
head(data_obs$Y_lag)
# check "data_obs"
head(data_obs) # starting 2003Q3 
tail(data_obs) # ending 2022Q3
#### OLS estimation ########################
OLS_out <- lm(MrP ~ Y + R + MrP_lag
              + dum_Q1 + dum_Q2 + dum_Q3,
              data = data_obs)
summary(OLS_out) 
names(OLS_out)
# coefficient of "Y"
beta_OLS <- OLS_out$coefficients[c(2, 3)]
print(beta_OLS)
#### TSLS estimation ########################
Y_est <- lm(Y ~ G + TAX + MrP_lag
            + CP_lag + IP_lag + R_lag
            + dum_Q1 + dum_Q2 + dum_Q3, 
            data = data_obs)
summary(Y_est) 
Y_hat <- Y_est$fitted.values
head(Y_hat)
print(c(length(Y), length(Y_hat)))
R_est <- lm(R ~ G + TAX + MrP_lag
            + CP_lag + IP_lag + R_lag
            + dum_Q1 + dum_Q2 + dum_Q3, 
            data = data_obs)
summary(R_est) 
R_hat <- R_est$fitted.values
head(R_hat)
data_TSLS <- data.frame(MrP = data_obs$MrP, 
                        Y_hat = Y_hat,
                        R_hat = R_hat,
                        MrP_lag = data_obs$MrP_lag,
                        dum_Q1 = data_obs$dum_Q1, 
                        dum_Q2 = data_obs$dum_Q2, 
                        dum_Q3 = data_obs$dum_Q3
                        )
TSLS_out <- lm(MrP ~ Y_hat + R_hat + MrP_lag
               + dum_Q1 + dum_Q2 + dum_Q3,  
               data = data_TSLS)
summary(TSLS_out) 
names(TSLS_out) 
print(TSLS_out$coefficients)
# coefficient of "Y_hat"
beta_TSLS <- TSLS_out$coefficients[c(2, 3)]
print(beta_TSLS)
print(beta_OLS)
#### LIML estimation ########################
# Some inputs of ivmodel() function are shown:
# Y: A numeric vector of outcomes.
# D: A vector of endogenous variables.
# Z: A matrix or data frame of instruments.
# X: A matrix or data frame of (exogenous) covariates.
#
# The model is: Y = X*alpha + D*beta + error
# E(error given X and Z) = 0
# Main interest is: beta.
# (1) dependent variable of "MrP" equation
Y_ivm = data_obs[, "MrP"]
# names(Y_ivm) <- "MrP" 
head(Y_ivm)
# (2) one endogenous explanatory variable in "MrP" equation,
D_ivm <- data_obs[, "Y"]
# names(D_ivm) <- "Y" 
head(D_ivm)
# (3) intruments for D_ivm = Y
# one instrument case
# Z_ivm <- data_obs[, "Y_lag"]
# names(Z_ivm) <- "Y_lag" 
# more than two instruments case
# The choice of "Z_ivm" is important for the output.
# Z_ivm <- data_obs[, c("Y_lag", "G", "TAX")]
# colnames(Z_ivm) <- c("Y_lag", "G", "TAX")
# Z_ivm <- data_obs[, c("CP", "IP", "G")]
# Z_ivm <- data_obs[, c("CP_lag", "IP_lag", "R_lag", "CP", "IP", "G", "TAX")]
# Z_ivm <- data_obs[, c("CP_lag", "IP_lag", "G", "TAX", "dum_Q1", "dum_Q2", "dum_Q3")]
# Z_ivm <- data_obs[, c("Y_lag", "R_lag", "G", "TAX")]
# do not use "G" for instruments due to use of "G <- Y - CP - IP".
Z_ivm <- data_obs[, c("Y_lag", "R_lag", "TAX")] 
head(Z_ivm)
# (4) exogenous explanatory variables in "MrP" equation,
# ivmodel::ivmodel() function has the limitation of dim(D) = 1.
# => "R" is treated as an exogenous explanatory variable.
# "R" is included in X_ivm. 
X_name <- c("R", "MrP_lag", "dum_Q1", "dum_Q2", "dum_Q3")
X_ivm <- data_obs[, X_name]
head(X_ivm)
#### input data to ivmodel() ###############################
# Y_ivm: "MrP"
# D_ivm: "Y"
head(cbind(Y_ivm, D_ivm, Z_ivm, X_ivm))
#### get "ivmodel" object ##################################
out_ivm <- ivmodel(Y = Y_ivm, D = D_ivm, Z = Z_ivm, X = X_ivm, 
                   intercept = TRUE,
                   beta0 = 0, alpha = 0.05,
                   k = c(0, 1),
                   manyweakSE = TRUE, heteroSE = FALSE, clusterID = NULL,
                   deltarange = NULL, na.action = na.omit
                   ) 
#### important notes #########################################
# (1) It is important to use "manyweakSE = TRUE".
# (2) The use of "k = c(0,3)" does not solve the problem.
# cheack the output of ivmodel()
print(out_ivm)
#### do LIML estimation #####################################
# The LIML estimator has a certain value of k.
out_LIML_ivm <- LIML(out_ivm,
                     beta0 = 0, alpha = 0.05,
                     manyweakSE = TRUE, heteroSE = FALSE, clusterID = NULL
                     )
# summary() does not work as expected.
summary(out_LIML_ivm) 
summary(out_LIML_ivm$point.est) 
print(out_LIML_ivm)
print(out_LIML_ivm$point.est)
print(out_LIML_ivm$point.est.other)
coef_out <- coef.ivmodel(out_LIML_ivm)
print(coef_out)
beta_LIML_ivm <- out_LIML_ivm$point.est
names(beta_LIML_ivm) <- "Y"
print(beta_LIML_ivm)
#### comparison of the coeffients of "Y" #####################
print(beta_OLS)
print(beta_TSLS)
print(beta_LIML_ivm[1])
#### Comments ################################################
# (1) ivmodel() function is very sensitive to the choice of instruments "Z".
# (2) It is difficult to find a good combination of instruments for the case of "MrP" equation.
# (3) ivmodel() function has the limitation of dim(D) = 1.
# (4) It is better to use TSLS estiomation for the case of dim(D) > 1. 
#### end #####################################################
q()
