#### scale check of y_graph
y_min <- min(y_graph)
y_max <- max(y_graph)
print(c(y_min, y_max))
#### plotting #####################################################
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 + 5), 
    ylim = c(80, 120), 
    xlab = "time (year)", 
    ylab = "food consumption per capita (dollars)", 
    main = "out-sample forecast, Kmenta Model, y = consump",
    sub = "in-sample = [1922 - 1941], out-sample = [1942 - 1945]")
## additional information to the graph
# abline(h = 0, col = "gray")
obs_last <- 1922 + 19  # year = 1941
text(obs_last, y_graph[(x_data == obs_last), 1], labels = "y_obs", pos = 4)
text(x_max, y_graph[(x_data == x_max), 2], labels = "y_fcst", pos = 4)
# pos = 1: below,  pos = 2: left,  pos = 3: up,  pos = 4: right.
#### save the graph into a file
# dev.copy(jpeg, file = "./outputs/forecast_out_sample_bimets_Kmenta_OK.jpg")
# dev.off()
#### end #############################################
#### R script ########################################
# Last modfied: 2024-12-09
# File name: forcast_out_sample_bimets_Kmenta_noGood.R
# Keywords: 
#   forcast_out_sample: out-sample forecast,
#   bimets: using "bimets" package,
#   Kments: the Kmenta model,  
#      using data set "Kmenta" in in "systemfit" package,
#   noGood: bimets does not work with the orinal Kmenta model,
#      using 3 equation of (Q_D, Q_S, Q_D = QS).
#            
# Written by: Hiroshi Murao
#
# Purpose: perfome the estimation of Kmenta model,
#      using bimets::ESTIMATION() function.  
#
# Data source:
#    Data set "Kmenta" in "systemfit" package.
# 
# List of variables:
#
# Q: food consumption per capita.
# P: ratio of food prices to general consumer prices.
# D: disposable income in constant dollars.
# F: ratio of preceding year's prices received 
#      by farmers to general consumer prices.
# A: time trend in years.
#
##########################################################
#### first things to do 
rm(list = ls())
setwd("C:/book_SIMU/data_works/Kmenta/R_prg")
library(systemfit) # Estimating Systems of Simultaneous Equations
library(bimets) # Time Series and Econometric Modeling
search()
#### get data
data("Kmenta")
head(Kmenta)
names(Kmenta)
#### make variables
consump <- Kmenta$consump       
price <- Kmenta$price     
income <- Kmenta$income     
farmPrice <- Kmenta$farmPrice 
trend <- Kmenta$trend
date_Y <- Kmenta$trend + 1921
consump_D <- Kmenta$consump
consump_S <- Kmenta$consump
#### convert data "Kmenta" for bimets ##########
# date_Y contains a numerical value such as 1922.
head(date_Y)
# convert 1922 to "1922-12-03"
# using as.character() function and as.Date() function
# convert, using as.character() function and as.Date() function
date_Ymd <- as.Date(as.character(date_Y), format = "%Y")
head(date_Ymd)
# date_Ymd contains a character value such as "1922-12-03".
data_obs <- as.data.frame(cbind(date_Ymd, Kmenta, consump_D, consump_S)) 
head(data_obs)   
#### convert "data_obs" to "myModelData" for bimets
myModelData <- lapply(as.list(xts(data_obs[, -1], 
               order.by = as.Date(data_obs[, 1],"%m/%d/%Y"))), 
               as.bimets)  
head(myModelData)  
# head(myModelData)
# tail(myModelData)
#### about Kmenta model ####################################
# A*y = B*x + C
# y: 2 endogenous variables => (consump, price)  
# x: 4 variables including "1", (income, farmPrice, trend, 1)     
#
# For bimets, two diffrent equations should have 
#   two different object names for "BEHAVIORAL". 
# We use the following variables: consump = consum_D = consump_S.
#### make myModel ####################################
myModel <- "
MODEL
COMMENT> Demand Equation
BEHAVIORAL> consump_D
EQ> consump_D = a12*price + b11*income + c1
COEFF> a12 b11 c1
COMMENT> Supply Equation 
BEHAVIORAL> consump_S
EQ> consump_S = a22*price + b22*farmPrice + b23*trend + c2 
COEFF> a22 b22 b23 c2 
COMMENT> Market Clearing 
IDENTITY> consump_D
EQ> consump_D = consump_S 
END
"
#### important notes
# (1) The use of the following commands makes "strange" outputs 
#     for forecast.
# 
# COMMENT> Market Clearing 
# IDENTITY> consump_D
# EQ> consump_D = consump_S 
#
# (2) We need the combination of (price, consump) for foreast. 
#
# COMMENT> Market Clearing 
# IDENTITY> consump_D
# EQ> consump_D = consump_S 
#### load "myModel" into "KmentaModel" #####################
KmentaModel <- LOAD_MODEL(modelText = myModel)
#### load "myModelData" into "KmentaModel" ################ 
KmentaModel <- LOAD_MODEL_DATA(KmentaModel, myModelData)
#### estimate the model ################################
# using ESTIMATE() function in "bimets" package 
IV_names <- c("1", "income", 
              "farmPrice", 
              "trend"
              )
obs_range <- c(1922, 1, 1941, 1)
est_range <- c(1922, 1, 1941, 1)
KmentaModel <- ESTIMATE(KmentaModel, 
           TSRANGE = est_range,           
           estTech = "IV", 
           IV = IV_names)
summary(KmentaModel) 
#### in-sample forecast ##############################################
fcst_range <- c(1923, 1, 1941, 1)
# 19 years from 1923 to 1941
# startining point = 1922 + 1
KmentaModel <- SIMULATE(KmentaModel,
              simType = 'FORECAST',
              TSRANGE = frcst_range,
              simConvergence = 0.01,
              simIterLimit = 100
              )
#### check the content of "KmentaModel"
names(KmentaModel)
names(KmentaModel$simulation)
str(KmentaModel$simulation)
#### get the output of "in-sample" forecast
GETDATE(KmentaModel$modelData$consump_D, format = "%Y-%m-%d")
obs_range <- c(1922, 1, 1941, 1)
est_range <- c(1922, 1, 1941, 1)
fcst_range <- c(1923, 1, 1941, 1)
# 19 years from 1923 to 1941
# The following line is important for plotting.
# "y" could be any endogenous variable.
y_fcst_part1 <- as.vector(fromBIMETStoTS(KmentaModel$simulation$consump_D))
print(y_fcst_part1)
print(length(y_fcst_part1))
TABIT(fromBIMETStoTS(KmentaModel$simulation$consump_D))
y_obs_raw <- as.vector(myModelData$consump_D)
head(y_obs_raw)
TABIT(myModelData$consump)
print(c(length(y_fcst_part1), length(y_obs_raw)))
# frcst_range = c(1923,1, 1941,1)
# obs_range <- c(1922, 1, 1941, 1)
# The following line is important for plotting.
y_obs_raw <- as.vector(myModelData$consump_D)
y_obs_part1 <- y_obs_raw[-c(1)] 
head(y_obs_raw)
head(y_obs_part1)
print(c(length(y_fcst_part1), length(y_obs_part1)))
y_part1 <- cbind(y_obs_part1, y_fcst_part1)
print(y_part1)
# There a different combination of data sets for "y_graph".
# Here we follow the bimets manual. 
#### make extention data for "out-sample" forecast #####
KmentaModel$modelData <- within(KmentaModel$modelData,{
  income = TSEXTEND(income, UPTO = c(1945,1), 
                     EXTMODE = "MYRATE", FACTOR = 1.02)
  farmPrice = TSEXTEND(farmPrice, UPTO = c(1945,1), 
                   EXTMODE = "MYRATE", FACTOR = 1.02)
  trend = TSEXTEND(trend, UPTO = c(1945,1), 
                   EXTMODE='LINEAR')
})
# FACTOR = 1.01 => 1% growth rate per period (one year for our case).
# check the content of extended data
TABIT(KmentaModel$modelData$income)
#### forecast of "out-sample" only ############################################
# Here we follow the bimets manual for merging two data sets.
fcst_range <- c(1942, 1, 1945, 1)
KmentaModel <- SIMULATE(KmentaModel,
              simType = 'FORECAST',
              TSRANGE = fcst_range,
              simConvergence = 0.00001,
              simIterLimit = 100
              )
#### check the result of simulation 
summary(KmentaModel)
names(KmentaModel)
names(Kmenta$simulation)
#### get the output of "out-sample" forecast
y_fcst_part2 <- as.vector(fromBIMETStoTS(KmentaModel$simulation$consump_D))
# 4 years from 1942 to 1945
print(y_fcst_part2)
#### preparations for plotting ############################################
# num_both = num_part1 + num_part2
num_part2 <- length(y_fcst_part2)
num_part1 <- length(y_fcst_part1)
print(c(num_part1, num_part2))
num_both <- num_part1 + num_part2
print(c(num_part1, num_part2, num_both))
y_obs_part2 <- matrix(NA, nrow = num_part2, 1)
y_part2 <- cbind(y_obs_part2, y_fcst_part2)
#### get y_graph, merging two data sets
y_graph <- rbind(y_part1, y_part2)
colnames(y_graph) <- c("y", "y_fcst")
print(y_graph)
#### get time horizon
num_H <- nrow(y_graph)
print(num_H)
x_data <- seq(1, num_H, by = 1)
x_data <- x_data + 1922
x_max <- max(x_data)
x_min <- min(x_data)
print(c(x_min, x_max))
# x_max = 2 times of nrow(y_graph) => Strange !!!!!!!!!!!!
#### scale ckeck of y_graph,
y_min <- min(y_graph)
y_max <- max(y_graph)
print(c(y_min, y_max))
#### plottinig ###################################################
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 + 10), 
    ylim = c(80, 120), 
    xlab = "time (year)", 
    ylab = "food consumption per capita (dollars)", 
    main = "out-sample forecast, Kmenta Model, consump",
    sub = "in-sample = [1922 - 1941], out-sample = [1942 - 1945]")
## additional information to the graph
# abline(h = 0, col = "gray")
obs_last <- 1922 + 19  # year = 1941
text(obs_last, y_graph[(x_data == obs_last), 1], labels = "y_obs", pos = 4)
text(x_max, y_graph[(x_data == x_max), 2], labels = "y_fcst", pos = 4)
# pos = 1: below,  pos = 2: left,  pos = 3: up,  pos = 4: right.
#### save the graph into a file
# dev.copy(jpeg, file = "./outputs/forecast_out_sample_bimets_Kmenta_noGood.jpg")
# dev.off()
#### Comments #######################################
# (1) The Kmenta model does not work well for forecast, 
#     for some software, including "bimets".
# (2) The use of the same Q for two endogenous variables 
#     becomes a problem to some or all kinds of software.
#### end #############################################
#### R script ########################################
# Last modfied: 2024-12-09
# File name: forcast_out_sample_bimets_Kmenta_noGood.R
# Keywords: 
#   forcast_out_sample: out-sample forecast,
#   bimets: using "bimets" package,
#   Kments: the Kmenta model,  
#      using data set "Kmenta" in in "systemfit" package,
#   noGood: bimets does not work with the orinal Kmenta model,
#      using 3 equation of (Q_D, Q_S, Q_D = QS).
#            
# Written by: Hiroshi Murao
#
# Purpose: perfome the estimation of Kmenta model,
#      using bimets::ESTIMATION() function.  
#
# Data source:
#    Data set "Kmenta" in "systemfit" package.
# 
# List of variables:
#
# Q: food consumption per capita.
# P: ratio of food prices to general consumer prices.
# D: disposable income in constant dollars.
# F: ratio of preceding year's prices received 
#      by farmers to general consumer prices.
# A: time trend in years.
#
##########################################################
#### first things to do 
rm(list = ls())
setwd("C:/book_SIMU/data_works/Kmenta/R_prg")
library(systemfit) # Estimating Systems of Simultaneous Equations
library(bimets) # Time Series and Econometric Modeling
search()
#### get data
data("Kmenta")
head(Kmenta)
names(Kmenta)
#### make variables
consump <- Kmenta$consump       
price <- Kmenta$price     
income <- Kmenta$income     
farmPrice <- Kmenta$farmPrice 
trend <- Kmenta$trend
date_Y <- Kmenta$trend + 1921
consump_D <- Kmenta$consump
consump_S <- Kmenta$consump
#### convert data "Kmenta" for bimets ##########
# date_Y contains a numerical value such as 1922.
head(date_Y)
# convert 1922 to "1922-12-03"
# using as.character() function and as.Date() function
# convert, using as.character() function and as.Date() function
date_Ymd <- as.Date(as.character(date_Y), format = "%Y")
head(date_Ymd)
# date_Ymd contains a character value such as "1922-12-03".
data_obs <- as.data.frame(cbind(date_Ymd, Kmenta, consump_D, consump_S)) 
head(data_obs)   
#### convert "data_obs" to "myModelData" for bimets
myModelData <- lapply(as.list(xts(data_obs[, -1], 
               order.by = as.Date(data_obs[, 1],"%m/%d/%Y"))), 
               as.bimets)  
head(myModelData)  
# head(myModelData)
# tail(myModelData)
#### about Kmenta model ####################################
# A*y = B*x + C
# y: 2 endogenous variables => (consump, price)  
# x: 4 variables including "1", (income, farmPrice, trend, 1)     
#
# For bimets, two diffrent equations should have 
#   two different object names for "BEHAVIORAL". 
# We use the following variables: consump = consum_D = consump_S.
#### make myModel ####################################
myModel <- "
MODEL
COMMENT> Demand Equation
BEHAVIORAL> consump_D
EQ> consump_D = a12*price + b11*income + c1
COEFF> a12 b11 c1
COMMENT> Supply Equation 
BEHAVIORAL> consump_S
EQ> consump_S = a22*price + b22*farmPrice + b23*trend + c2 
COEFF> a22 b22 b23 c2 
COMMENT> Market Clearing 
IDENTITY> consump_D
EQ> consump_D = consump_S 
END
"
#### important notes
# (1) The use of the following commands makes "strange" outputs 
#     for forecast.
# 
# COMMENT> Market Clearing 
# IDENTITY> consump_D
# EQ> consump_D = consump_S 
#
# (2) We need the combination of (price, consump) for foreast. 
#
# COMMENT> Market Clearing 
# IDENTITY> consump_D
# EQ> consump_D = consump_S 
#### load "myModel" into "KmentaModel" #####################
KmentaModel <- LOAD_MODEL(modelText = myModel)
#### load "myModelData" into "KmentaModel" ################ 
KmentaModel <- LOAD_MODEL_DATA(KmentaModel, myModelData)
#### estimate the model ################################
# using ESTIMATE() function in "bimets" package 
IV_names <- c("1", "income", 
              "farmPrice", 
              "trend"
              )
obs_range <- c(1922, 1, 1941, 1)
est_range <- c(1922, 1, 1941, 1)
KmentaModel <- ESTIMATE(KmentaModel, 
           TSRANGE = est_range,           
           estTech = "IV", 
           IV = IV_names)
summary(KmentaModel) 
#### in-sample forecast ##############################################
fcst_range <- c(1923, 1, 1941, 1)
# 19 years from 1923 to 1941
# startining point = 1922 + 1
KmentaModel <- SIMULATE(KmentaModel,
              simType = 'FORECAST',
              TSRANGE = fcst_range,
              simConvergence = 0.01,
              simIterLimit = 100
              )
#### check the content of "KmentaModel"
names(KmentaModel)
names(KmentaModel$simulation)
str(KmentaModel$simulation)
#### get the output of "in-sample" forecast
GETDATE(KmentaModel$modelData$consump_D, format = "%Y-%m-%d")
obs_range <- c(1922, 1, 1941, 1)
est_range <- c(1922, 1, 1941, 1)
fcst_range <- c(1923, 1, 1941, 1)
# 19 years from 1923 to 1941
# The following line is important for plotting.
# "y" could be any endogenous variable.
y_fcst_part1 <- as.vector(fromBIMETStoTS(KmentaModel$simulation$consump_D))
print(y_fcst_part1)
print(length(y_fcst_part1))
TABIT(fromBIMETStoTS(KmentaModel$simulation$consump_D))
y_obs_raw <- as.vector(myModelData$consump_D)
head(y_obs_raw)
TABIT(myModelData$consump)
print(c(length(y_fcst_part1), length(y_obs_raw)))
# frcst_range = c(1923,1, 1941,1)
# obs_range <- c(1922, 1, 1941, 1)
# The following line is important for plotting.
y_obs_raw <- as.vector(myModelData$consump_D)
y_obs_part1 <- y_obs_raw[-c(1)] 
head(y_obs_raw)
head(y_obs_part1)
print(c(length(y_fcst_part1), length(y_obs_part1)))
y_part1 <- cbind(y_obs_part1, y_fcst_part1)
print(y_part1)
# There a different combination of data sets for "y_graph".
# Here we follow the bimets manual. 
#### make extention data for "out-sample" forecast #####
KmentaModel$modelData <- within(KmentaModel$modelData,{
  income = TSEXTEND(income, UPTO = c(1945,1), 
                     EXTMODE = "MYRATE", FACTOR = 1.02)
  farmPrice = TSEXTEND(farmPrice, UPTO = c(1945,1), 
                   EXTMODE = "MYRATE", FACTOR = 1.02)
  trend = TSEXTEND(trend, UPTO = c(1945,1), 
                   EXTMODE='LINEAR')
})
# FACTOR = 1.01 => 1% growth rate per period (one year for our case).
# check the content of extended data
TABIT(KmentaModel$modelData$income)
#### forecast of "out-sample" only ############################################
# Here we follow the bimets manual for merging two data sets.
fcst_range <- c(1942, 1, 1945, 1)
KmentaModel <- SIMULATE(KmentaModel,
              simType = 'FORECAST',
              TSRANGE = fcst_range,
              simConvergence = 0.00001,
              simIterLimit = 100
              )
#### check the result of simulation 
summary(KmentaModel)
names(KmentaModel)
names(Kmenta$simulation)
#### get the output of "out-sample" forecast
y_fcst_part2 <- as.vector(fromBIMETStoTS(KmentaModel$simulation$consump_D))
# 4 years from 1942 to 1945
print(y_fcst_part2)
#### preparations for plotting ############################################
# num_both = num_part1 + num_part2
num_part2 <- length(y_fcst_part2)
num_part1 <- length(y_fcst_part1)
print(c(num_part1, num_part2))
num_both <- num_part1 + num_part2
print(c(num_part1, num_part2, num_both))
y_obs_part2 <- matrix(NA, nrow = num_part2, 1)
y_part2 <- cbind(y_obs_part2, y_fcst_part2)
#### get y_graph, merging two data sets
y_graph <- rbind(y_part1, y_part2)
colnames(y_graph) <- c("y", "y_fcst")
print(y_graph)
#### get time horizon
num_H <- nrow(y_graph)
print(num_H)
x_data <- seq(1, num_H, by = 1)
x_data <- x_data + 1922
x_max <- max(x_data)
x_min <- min(x_data)
print(c(x_min, x_max))
# x_max = 2 times of nrow(y_graph) => Strange !!!!!!!!!!!!
#### scale ckeck of y_graph,
y_min <- min(y_graph)
y_max <- max(y_graph)
print(c(y_min, y_max))
#### plottinig ###################################################
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 + 10), 
    ylim = c(80, 120), 
    xlab = "time (year)", 
    ylab = "food consumption per capita (dollars)", 
    main = "out-sample forecast, Kmenta Model, consump",
    sub = "in-sample = [1922 - 1941], out-sample = [1942 - 1945]")
## additional information to the graph
# abline(h = 0, col = "gray")
obs_last <- 1922 + 19  # year = 1941
text(obs_last, y_graph[(x_data == obs_last), 1], labels = "y_obs", pos = 4)
text(x_max, y_graph[(x_data == x_max), 2], labels = "y_fcst", pos = 4)
# pos = 1: below,  pos = 2: left,  pos = 3: up,  pos = 4: right.
#### save the graph into a file
# dev.copy(jpeg, file = "./outputs/forecast_out_sample_bimets_Kmenta_noGood.jpg")
# dev.off()
#### Comments #######################################
# (1) The Kmenta model does not work well for forecast, 
#     for some software, including "bimets".
# (2) The use of the same Q for two endogenous variables 
#     becomes a problem to some or all kinds of software.
#### end #############################################
q()
