#   )
print(plot)
ggsave("figures/fig_fitforc_10y.pdf", plot = plot,width = 5.72, height = 3.50)
plotforc3m <- data.frame(data.frame(ACM$Date,ACM$n_1,ACM$n_1for,
ACM$n_1fit
) )
# Convert the data frame to a data.table
setDT(plotforc3m)
# Create the n_40plot series
plotforc3m[, ACM.n_1plot := ifelse(ACM.Date <= as.Date("2005-12-01"), ACM.n_1fit, ACM.n_1for)]
plot<-ggplot(plotforc3m, aes(x = ACM.Date)) +
geom_line(aes(y = ACM.n_1, color = "n_1"), size = 1) +
geom_line(aes(y = ACM.n_1plot, color = "Standard Model"), linetype = "dashed", size = 1) +
labs(x = "Date", y = "Fit and Forecast", title = "3-month Yields") +
scale_color_manual(values = c("n_1" = "blue", "Standard Model" = "red")) +
theme_minimal() +
theme(
legend.position = "top",
plot.title = element_text(hjust = 0.5)  # Center the title
) +
guides(color = guide_legend(title = ""))
plot <- plot +
geom_vline(xintercept = as.Date("2005-12-01"), color = "black")
# Shading the area from 2006 onwards
# plot <- plot +
#   geom_rect(
#     xmin = as.Date("2006-01-01"),
#     xmax = max(plotforc3m$ACM.Date),
#     ymin = -Inf,
#     ymax = Inf,
#     fill = "gray",
#     alpha = 0.01,
#     color = NA  # Set the outline color to match the fill color
#   )
print(plot)
ggsave("figures/fig_fitforc_3m.pdf", plot = plot,width = 5.72, height = 3.50)
################################################################################
###
### AUTHOR:       Carlo Favero
### DATE:         Fall 2023
### DESCRIPTION:  This script creates a database with all TS data
### OUTPUT:       databases for the term  for the term
###
###################################################################################
## -----------------------------------------------------------------------------
rm(list=ls())                                                                   # Clear the environment
## -----------------------------------------------------------------------------
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))                     # Set the appropriate environment to work in every machine
## -----------------------------------------------------------------------------
listofpackages <- c("tidyverse","ellipse","reshape2","xts","xlsx","readxl",     # Install and load the required packages
"quantmod")
source("../2024_ex1/preparePackages.R")
preparePackages(listofpackages)
## -----------------------------------------------------------------------------
# Nelson, Svensson, Siegel yield curve parametrization
nss <- function(n, params) {
Beta0 <- params$BETA0
Beta1 <- params$BETA1
Beta2 <- params$BETA2
Beta3 <- params$BETA3
Tau1  <- params$TAU1
Tau2  <- params$TAU2
Tau1
return(Beta0 + Beta1 * (1 - exp(-n / Tau1)) / (n / Tau1) +
Beta2 * ((1 - exp(-n / Tau1)) / (n / Tau1) - exp(-n / Tau1)) +
Beta3 * ((1 - exp(-n / Tau2)) / (n / Tau2) - exp(-n / Tau2)))
}
### YIELDS
### long-term yields from Gurkanyak
### use the betas and the taus in the file to construct yields at all needed maturities
### name them in terms of quarters and add them to the database
### the file with all the info is yields_GURK.xlsx
yields<- read_excel('data/yields_GURK.xlsx', sheet=1, col_names=TRUE)
dates <-seq(as.Date("1971-12-01"),length=nrow(yields), by="quarters")
yields.ts <- xts(yields, order.by=dates)
yields.ts <- yields.ts[, -1]
# Take just the parameters
params.cols <- grep("^(BETA|TAU)", colnames(yields.ts))
params.ts <- yields.ts[, params.cols]
save(params.ts,file="data/yields_GURK_params.rda")
write.csv(data.frame(params.ts), "data/yields_GURK_params.csv", row.names = TRUE)
# Take the yields and save their names as "nX" where X is a number
yields.cols <- grep("^SVENY", colnames(yields.ts))
yields.ts <- yields.ts[, yields.cols]
colnames(yields.ts) <- gsub("^SVENY", "y", colnames(yields.ts))
new_colnames <- sapply(str_extract_all(colnames(yields.ts), "\\d+"), function(x) {
paste0("y", as.integer(x) * 4)
})
colnames(yields.ts) <- new_colnames
# Since computations are a nightmare with time series, we convert everything to numeric
# and do everything with this class
Beta0 <- as.numeric(params.ts$BETA0)
Beta1 <- as.numeric(params.ts$BETA1)
Beta2 <- as.numeric(params.ts$BETA2)
Beta3 <- as.numeric(params.ts$BETA3)
Tau1  <- as.numeric(params.ts$TAU1)
Tau2  <- as.numeric(params.ts$TAU2)
N=60
ttm = c(1:N)*(1/4)
maturities<-c(1:N)
for (n in maturities){
new_yield <-(Beta0 + Beta1 * (1 - exp(-0.25*n / Tau1)) / (0.25*n / Tau1) +
Beta2 * ((1 - exp(-0.25*n / Tau1)) / (0.25*n / Tau1) - exp(-0.25*n / Tau1)) +
Beta3 * ((1 - exp(-0.25*n / Tau2)) / (0.25*n / Tau2) - exp(-0.25*n / Tau2)))
# Again to yields
new_yield <- xts(new_yield, order.by=index(yields.ts))
# Fancy header name
colname   <- paste0("yNS", n)
# Merge it with respect to the old dataset
yields.ts <- merge(yields.ts, new_yield, join="left")
colnames(yields.ts)[ncol(yields.ts)] <- colname
}
nss_yields<-0.01*yields.ts[,16:75]
colnames(nss_yields) <- paste("n",seq(1,N),sep = "_")
rm(list="params.ts","yields","yields.ts","new_yield")
save(nss_yields,file="data/yields1.rda")
## -----------------------------------------------------------------------------
rm(list=ls())                                                                   # Clear the environment
## -----------------------------------------------------------------------------
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))                     # Set the appropriate environment to work                                                                                # in every machine
## -----------------------------------------------------------------------------
listofpackages <- c("psych","pracma","zoo","lubridate","dplyr",
"matrixcalc","matlib","xlsx","readxl","readr","xts","tidyverse","vars","systemfit","urca")
source("../2024_ex1/preparePackages.R")
preparePackages(listofpackages)
rm(listofpackages, preparePackages)
## -----------------------------------------------------------------------------
## Load NS data
load('data/yields1.rda')
#identify the maturity, in this case 60 quarters one period is one quarter.
N=60
ttm = c(1:N)*(1/4)
maturities<-c(1:N)
nss_yields<-nss_yields["1980-03-01/2023-06-01"]
dates <- index(nss_yields)
dates_lag <-seq(as.Date("1980-06-01"),length=nrow(nss_yields)-1, by="quarters")
########################
# K Principal Components
########################
K=5
demeanedyield = scale(nss_yields, scale = FALSE)
k = prcomp((demeanedyield))
yieldPCs = demeanedyield %*% (k$rotation)
# get the principal components in .xts format and add them to yields.ts
X<- xts(yieldPCs[,1:K], order.by=dates)
X.df <- data.frame(
Date = index(X),X)
save(X.df,file="data/output/X_ACM.rda")
nss_yields <- merge(nss_yields, X, join = 'outer')
################################################################
#Step 1 - estimate VAR(1) for the time series of pricing factors
################################################################
#using.xts, this needs the package "vars"
#my_var_model <- VAR(X, p = 1, type = "none")
my_var_model <- VAR(X, p = 1, type = "const")
#cointegration analysis
summary(my_var_model)
result <- ca.jo(X, type = "eigen", ecdet = "const",K=2)
summary(result)
coeff_matrices <- coef(my_var_model)
phi <- rbind(coeff_matrices$PC1[1:5,1],coeff_matrices$PC2[1:5,1],
coeff_matrices$PC3[1:5,1],coeff_matrices$PC4[1:5,1],
coeff_matrices$PC5[1:5,1])
mu<- rbind(coeff_matrices$PC1[6,1],coeff_matrices$PC2[6,1],
coeff_matrices$PC3[6,1],coeff_matrices$PC4[6,1],
coeff_matrices$PC5[6,1])
resmat <- residuals(my_var_model)
Sigma <-  t(resmat)%*%resmat/nrow(resmat)
#merge residuals into .xts database
colnames(resmat) <- paste("v",colnames(resmat),sep = "_")
resmat.ts <- xts(resmat, order.by=dates_lag)
nss_yields <- merge(nss_yields, resmat.ts, join = 'outer')
####################################
#Step 2 - regress log excess returns
####################################
#generate excess returns
#generate hpr and excess hpr
nss_yields$hpr4 <- (1/4)*lag.xts(nss_yields$n_4) - (4/4-1/4)*(nss_yields$n_3-lag.xts(nss_yields$n_4))
nss_yields$hpr8 <- (1/4)*lag.xts(nss_yields$n_8) - (8/4-1/4)*(nss_yields$n_7-lag.xts(nss_yields$n_8))
nss_yields$hpr12<- (1/4)*lag.xts(nss_yields$n_12)- (12/4-1/4)*(nss_yields$n_11-lag.xts(nss_yields$n_12))
nss_yields$hpr16<- (1/4)*lag.xts(nss_yields$n_16)- (16/4-1/4)*(nss_yields$n_15-lag.xts(nss_yields$n_16))
nss_yields$hpr20<- (1/4)*lag.xts(nss_yields$n_20)- (20/4-1/4)*(nss_yields$n_19-lag.xts(nss_yields$n_20))
nss_yields$hpr28<- (1/4)*lag.xts(nss_yields$n_28)- (28/4-1/4)*(nss_yields$n_27-lag.xts(nss_yields$n_28))
nss_yields$hpr40<- (1/4)*lag.xts(nss_yields$n_40)- (40/4-1/4)*(nss_yields$n_39-lag.xts(nss_yields$n_40))
nss_yields$hpr60<- (1/4)*lag.xts(nss_yields$n_60)- (60/4-1/4)*(nss_yields$n_59-lag.xts(nss_yields$n_60))
nss_yields$hper4=nss_yields$hpr4-(1/4)*lag.xts(nss_yields$n_1)
nss_yields$hper8=nss_yields$hpr8-(1/4)*lag.xts(nss_yields$n_1)
nss_yields$hper12=nss_yields$hpr12-(1/4)*lag.xts(nss_yields$n_1)
nss_yields$hper16=nss_yields$hpr16-(1/4)*lag.xts(nss_yields$n_1)
nss_yields$hper20=nss_yields$hpr20-(1/4)*lag.xts(nss_yields$n_1)
nss_yields$hper28=nss_yields$hpr28-(1/12)*lag.xts(nss_yields$n_1)
nss_yields$hper40=nss_yields$hpr40-(1/4)*lag.xts(nss_yields$n_1)
nss_yields$hper60=nss_yields$hpr60-(1/4)*lag.xts(nss_yields$n_1)
#lags of factors
nss_yields$PC1_lag=lag.xts(nss_yields$PC1)
nss_yields$PC2_lag=lag.xts(nss_yields$PC2)
nss_yields$PC3_lag=lag.xts(nss_yields$PC3)
nss_yields$PC4_lag=lag.xts(nss_yields$PC4)
nss_yields$PC5_lag=lag.xts(nss_yields$PC5)
###----------------------
###System estimation for excess returns
###----------------------
eqhper_4r <- hper4 ~  PC1_lag+ PC2_lag+PC3_lag+PC4_lag+PC5_lag
eqhper_8r <- hper8 ~  PC1_lag+ PC2_lag+PC3_lag+PC4_lag+PC5_lag
eqhper_12r <- hper12 ~  PC1_lag+ PC2_lag+PC3_lag+PC4_lag+PC5_lag
eqhper_16r <- hper16 ~  PC1_lag+ PC2_lag+PC3_lag+PC4_lag+PC5_lag
eqhper_20r <- hper20 ~  PC1_lag+ PC2_lag+PC3_lag+PC4_lag+PC5_lag
eqhper_28r <- hper28 ~  PC1_lag+ PC2_lag+PC3_lag+PC4_lag+PC5_lag
eqhper_40r <- hper40 ~  PC1_lag+ PC2_lag+PC3_lag+PC4_lag+PC5_lag
eqhper_60r <- hper60 ~  PC1_lag+ PC2_lag+PC3_lag+PC4_lag+PC5_lag
system_hper_r <- list(eqhper_4r,eqhper_8r,eqhper_12r,eqhper_16r,eqhper_20r,eqhper_28r,
eqhper_40r,eqhper_60r)
eqhper_4 <- hper4 ~  v_PC1+ v_PC2+v_PC3+v_PC4+v_PC5+
+PC1_lag+ PC2_lag+PC3_lag+PC4_lag+PC5_lag
eqhper_8 <- hper8 ~  v_PC1+ v_PC2+v_PC3+v_PC4+v_PC5+
+PC1_lag+ PC2_lag+PC3_lag+PC4_lag+PC5_lag
eqhper_12 <- hper12 ~  v_PC1+ v_PC2+v_PC3+v_PC4+v_PC5+
+PC1_lag+ PC2_lag+PC3_lag+PC4_lag+PC5_lag
eqhper_16 <- hper16 ~  v_PC1+ v_PC2+v_PC3+v_PC4+v_PC5+
+PC1_lag+ PC2_lag+PC3_lag+PC4_lag+PC5_lag
eqhper_20 <- hper20 ~  v_PC1+ v_PC2+v_PC3+v_PC4+v_PC5+
+PC1_lag+ PC2_lag+PC3_lag+PC4_lag+PC5_lag
eqhper_28 <- hper28 ~  v_PC1+ v_PC2+v_PC3+v_PC4+v_PC5+
+PC1_lag+ PC2_lag+PC3_lag+PC4_lag+PC5_lag
eqhper_40 <- hper40 ~  v_PC1+ v_PC2+v_PC3+v_PC4+v_PC5+
+PC1_lag+ PC2_lag+PC3_lag+PC4_lag+PC5_lag
eqhper_60 <- hper60 ~ v_PC1+ v_PC2+v_PC3+v_PC4+v_PC5+
+PC1_lag+ PC2_lag+PC3_lag+PC4_lag+PC5_lag
system_hper_unr <- list(eqhper_4,eqhper_8,eqhper_12,eqhper_16,eqhper_20,eqhper_28,
eqhper_40,eqhper_60)
#estimation
system_hper_r_fit <- systemfit( system_hper_r, "OLS", data = nss_yields )
summary(system_hper_r_fit)
system_hper_unr_fit <- systemfit( system_hper_unr, "OLS", data = nss_yields )
summary(system_hper_unr_fit)
coeffmat_sys <- coef(system_hper_unr_fit)
res_sys<- residuals(system_hper_unr_fit)
res_sys<-res_sys[-1,] #first row contains NAs because of the lag
E<-as.matrix(res_sys)
sigmasq_ret = tr(t(E) %*% E)/(nrow(E)*ncol(E))
#Sigma_v1<-  (t(E1)%*%E1)/(nrow(E1))
#sigmasq_ret1 <- as.matrix(diag(Sigma_v1))
m1=matrix(coeffmat_sys,nrow=8,ncol=11,byrow=TRUE)
#m1 <- matrix(coeff_matsys, nrow = 11, ncol = 8)
a <- m1[,1]
beta <- t(m1[,2:(K+1)])
c <- t(m1[,(K+2):ncol(m1)])
#Step 3 - Get lambdas and Bond prices
quad_form <- function(b){
c(outer(b,b))
}
BStar <- t(apply(t(beta),1,quad_form))
lambda0 <- Ginv(t(beta)) %*% (a + 0.5*((BStar %*% c(Sigma)) + sigmasq_ret))
lambda1 <- Ginv(t(beta)) %*% t(c)
#Recursions for bond pricing
A <- rep(0,N)
B <- matrix(0,K,N)
nss_yields$rf<-0.25*nss_yields$n_1
mod3 <- lm(nss_yields$rf ~ X)
delta0 <- coef(mod3)[1]
delta1 <- coef(mod3)[-1]
A[1] = -delta0
B[, 1] = - delta1
for (i in 2:N){
A[i]  = A[i-1] + t(B[, i-1]) %*% (mu-lambda0) + 0.5 * (t(B[, i-1]) %*% Sigma %*% B[, i-1] + sigmasq_ret) - delta0
B[,i] = t(B[, i-1]) %*% (phi - lambda1) - delta1
}
#Recursions for Risk Free Yields
A_rf <- rep(0,N)
B_rf <- matrix(0,K,N)
A_rf[1] = -delta0
B_rf[, 1] = - delta1
lambda0 = matrix(0,K,1)
lambda1 = matrix(0,K,K)
for (i in 2:N){
A_rf[i]  = A_rf[i-1] + t(B_rf[, i-1]) %*% (mu-lambda0) + 0.5 * (t(B_rf[, i-1]) %*% Sigma %*% B_rf[, i-1] + sigmasq_ret) - delta0
B_rf[,i] = t(B_rf[, i-1]) %*% (phi - lambda1) - delta1
}
#FITTED YIELDS
fittedLogPrices = t(as.vector(t(A)) + t(B) %*% t(X))
fittedYields = - t(t(fittedLogPrices) / ttm)
fittedYields.ts <- xts(fittedYields, order.by=dates)
#save data for plots
plotactfit10Y <- data.frame(
Date = index(nss_yields),
Actual = nss_yields[, 40],
Fitted = fittedYields[, 40]
)
plotactfit3m <- data.frame(
Date = index(nss_yields),
Actual = nss_yields[, 1],
Fitted = fittedYields[, 1]
)
#RISK NEUTRAL YIELDS
RiskFreeLogPrices = t(as.vector(t(A_rf)) + t(B_rf) %*% t(X))
RiskFreeYields = - t(t(RiskFreeLogPrices) / ttm)
#TERM PREMIA
termpremia = fittedYields - RiskFreeYields
plot_tp <- data.frame(
Date = index(nss_yields),
TP40 = termpremia[, 40],
TP4 = termpremia[, 4]
)
#save(nss_yields, file="data/cf_nss_yields.rda")
###----------------------
###Forecasting
###----------------------
X_F<-X["1980-03-01/2005-12-01"]
var_forc <- VAR(X_F, p = 1, type = "const")
summary(var_forc)
X_FOR <- predict(var_forc, n.ahead = 70)
X_PRED<-cbind(X_FOR$fcst$PC1[,1],X_FOR$fcst$PC2[,1],X_FOR$fcst$PC3[,1],
X_FOR$fcst$PC4[,1],X_FOR$fcst$PC5[,1])
n_1_pred<- -(A[1]+t(B[, 1])%*%t(X_PRED))/ttm[1]
n_40_pred<- -(A[40]+t(B[, 40])%*%t(X_PRED))/ttm[40]
dates_for <-seq(as.Date("2006-03-01"),length=70, by="quarters")
n_1_pred.ts<- xts(t(n_1_pred), order.by=dates_for)
n_40_pred.ts<- xts(t(n_40_pred), order.by=dates_for)
n_1_pred.ts <- merge(n_1_pred.ts, nss_yields$n_1, join = 'outer')
n_40_pred.ts <- merge(n_40_pred.ts, nss_yields$n_40, join = 'outer')
colnames(n_40_pred.ts) <- c("n_40for","n_40")
colnames(n_1_pred.ts) <- c("n_1for","n_1")
plotforc3m <- data.frame(
Date = index(n_1_pred.ts),
Actual = coredata(n_1_pred.ts$n_1),
Forecast =coredata(n_1_pred.ts$n_1for) )
plotforc10y <- data.frame(
Date = index(n_40_pred.ts),
Actual = coredata(n_40_pred.ts$n_40),
Forecast =coredata(n_40_pred.ts$n_40for) )
# saving the data from the standard model
ACM <- data.frame(plotforc10y,plotforc3m[,2:3],
plot_tp[,2:3],plotactfit10Y[,3],plotactfit3m[,3]
)
colnames(ACM) <- c("Date","n_40","n_40for","n_1","n_1for","TP40","TP4","n_40fit","n_1fit")
save(ACM,file="data/output/ACM.rda")
## -----------------------------------------------------------------------------
rm(list=ls())                                                                   # Clear the environment
## -----------------------------------------------------------------------------
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))                     # Set the appropriate environment to work                                                                                # in every machine
## -----------------------------------------------------------------------------
listofpackages <- c("psych","pracma","zoo","lubridate","dplyr",
"matrixcalc","matlib","xlsx","readxl","readr","xts","tidyverse","vars","systemfit")
source("../2024_ex1/preparePackages.R")
preparePackages(listofpackages)
rm(listofpackages, preparePackages)
## -----------------------------------------------------------------------------
## Load NS data
load('data/yields1.rda')
#identify the maturity, in this case 60 quarters one period is one quarter.
N=60
ttm = c(1:N)*(1/4)
maturities<-c(1:N)
nss_yields<-nss_yields["1980-03-01/2023-06-01"]
dates <- index(nss_yields)
dates_lag <-seq(as.Date("1980-06-01"),length=nrow(nss_yields)-1, by="quarters")
#########################################
# Figure 1-2 yields and returns
#########################################
ytm_df <- data.frame(dates = index(nss_yields), Data = coredata(nss_yields))
col_names <- c("date", colnames(nss_yields))
colnames(ytm_df) <- col_names
Yields <- as.matrix(ytm_df[,2:61])
T = nrow(Yields)
logPrices = t(t(-Yields)*ttm)
rf =  -1*logPrices[1:T,1]
rx =  logPrices[2:T, 1:(N-1)] - logPrices[1:(T-1), 2:N] - rf[1:T-1]
rx_df<- data.frame(dates = dates_lag, Data = coredata(rx))
colnames_er<- paste("er",seq(2,N),sep = "_")
colnames(rx_df)<- c("date", colnames_er)
long_data <- ytm_df %>%
gather(key = "Variable", value = "Value", starts_with("n_"))
# Set the levels for the "Variable" factor to ensure the desired color mapping
long_data$Variable <- factor(long_data$Variable, levels = unique(long_data$Variable))
# Define a custom color palette from pale blue to dark blue for the y variables
color_palette <- colorRampPalette(c("lightblue", "darkblue"))(length(unique(long_data$Variable)))
plot <- ggplot(long_data, aes(x = date, y = Value, color = Variable)) +
geom_line(size = 1) +
labs(title = "Common Trends in Yields",
x = "Time", y = "Value") +
scale_color_manual(values = setNames(color_palette, levels(long_data$Variable))) +  # Manually set colors
#geom_line(data = subset_df, aes(x = Dates, y = y3m_cc), color = "red", size = 1) +  # Plot x in red
theme_minimal() +
theme(axis.line = element_line(color = "black")) +
scale_x_date(date_breaks = "5 year", date_labels = "%Y") +
theme(legend.text = element_text(size = 6),  # Adjust the font size in the legend
legend.key.size = unit(0.1, "cm"))  # Specify a smaller legend key size (1 cm)
print(plot)
ggsave(filename = "figures/Figure_1.pdf", plot = plot, device = "pdf",width = 5.72,height=3.12)
rm(long_data)
long_data <- rx_df %>%
gather(key = "Variable", value = "Value", starts_with("er_"))
# Set the levels for the "Variable" factor to ensure the desired color mapping
long_data$Variable <- factor(long_data$Variable, levels = unique(long_data$Variable))
# Define a custom color palette from pale blue to dark blue for the y variables
color_palette <- colorRampPalette(c("lightblue", "darkblue"))(length(unique(long_data$Variable)))
plot <- ggplot(long_data, aes(x = date, y = Value, color = Variable)) +
geom_line(size = 1) +
labs(title = "Common Cycles in Excess Period Holding Returns",
x = "Time", y = "Value") +
scale_color_manual(values = setNames(color_palette, levels(long_data$Variable))) +  # Manually set colors
#geom_line(data = subset_df, aes(x = Dates, y = y3m_cc), color = "red", size = 1) +  # Plot x in red
theme_minimal() +
theme(axis.line = element_line(color = "black")) +
scale_x_date(date_breaks = "5 year", date_labels = "%Y") +
theme(legend.text = element_text(size = 6),  # Adjust the font size in the legend
legend.key.size = unit(0.1, "cm"))  # Specify a smaller legend key size (1 cm)
print(plot)
ggsave(filename = "figures/Figure_2.pdf", plot = plot, device = "pdf",width = 5.72,height=3.12)
rm(plot,rx,rx_df,ytm_df,long_data,logPrices,Yields)
#########################################
# Factors
#########################################
load('data/output/X_ACM.rda')
load('data/output/ACM.rda')
X.df_long <- gather(X.df, key = "Variable", value = "Value", -Date)
# Create the time-series plot with a legend
plot <-ggplot(data = X.df_long, aes(x = Date, y = Value, color = Variable)) +
geom_line() +
labs(title = "The first five PC extracted from yields",
x = "Date",
y = "Value") +
scale_color_manual(values = c("PC1" = "red", "PC2" = "blue", "PC3" = "green", "PC4" = "purple", "PC5" = "orange")) +
theme_minimal() +
theme(legend.title = element_blank())
print(plot)
ggsave("figures/fig_X.pdf", plot = plot,width = 5.72, height = 3.12)
#######################
# TERM PREMIA
#######################
plot<- ggplot(ACM, aes(x = Date)) +
geom_line(aes(y = TP40, color = "TP40"), size = 1) +
geom_line(aes(y = TP4, color = "TP4"), linetype = "dashed", size = 1) +
labs(x = "Date", y = "Term PRemia", title = "1Y and 10Y Term Premia in the standard model") +
scale_color_manual(values = c("TP40" = "blue", "TP4" = "red")) +
theme_minimal() +
theme(
legend.position = "top",
plot.title = element_text(hjust = 0.5)  # Center the title
) +
guides(color = guide_legend(title = "Series"))
print(plot)
ggsave("figures/fig_ACM_tp.pdf", plot = plot)
#######################
# FIT and FORECAST
#######################
plotforc10y <- data.frame(data.frame(ACM$Date,ACM$n_40,ACM$n_40for,ACM$n_40fit
) )
library(data.table)
# Convert the data frame to a data.table
setDT(plotforc10y)
# Create the n_40plot series
plotforc10y[, ACM.n_40plot := ifelse(ACM.Date <= as.Date("2005-12-01"), ACM.n_40fit, ACM.n_40for)]
plot<-ggplot(plotforc10y, aes(x = ACM.Date)) +
geom_line(aes(y = ACM.n_40, color = "n_40"), size = 1) +
geom_line(aes(y = ACM.n_40plot, color = "Standard Model"), linetype = "dashed", size = 1) +
labs(x = "Date", y = "Fit and Forecast", title = "10-year Yields") +
scale_color_manual(values = c("n_40" = "blue", "Standard Model" = "red")) +
theme_minimal() +
theme(
legend.position = "top",
plot.title = element_text(hjust = 0.5)  # Center the title
) +
guides(color = guide_legend(title = ""))
plot <- plot +
geom_vline(xintercept = as.Date("2005-12-01"), color = "black")
# Shading the area from 2006 onwards
# plot <- plot +
#   geom_rect(
#     xmin = as.Date("2006-01-01"),
#     xmax = max(plotforc10y$ACM.Date),
#     ymin = -Inf,
#     ymax = Inf,
#     fill = "gray",
#     alpha = 0.01,
#     color = NA  # Set the outline color to match the fill color
#   )
print(plot)
ggsave("figures/fig_fitforc_10y.pdf", plot = plot,width = 5.72, height = 3.50)
plotforc3m <- data.frame(data.frame(ACM$Date,ACM$n_1,ACM$n_1for,
ACM$n_1fit
) )
# Convert the data frame to a data.table
setDT(plotforc3m)
# Create the n_40plot series
plotforc3m[, ACM.n_1plot := ifelse(ACM.Date <= as.Date("2005-12-01"), ACM.n_1fit, ACM.n_1for)]
plot<-ggplot(plotforc3m, aes(x = ACM.Date)) +
geom_line(aes(y = ACM.n_1, color = "n_1"), size = 1) +
geom_line(aes(y = ACM.n_1plot, color = "Standard Model"), linetype = "dashed", size = 1) +
labs(x = "Date", y = "Fit and Forecast", title = "3-month Yields") +
scale_color_manual(values = c("n_1" = "blue", "Standard Model" = "red")) +
theme_minimal() +
theme(
legend.position = "top",
plot.title = element_text(hjust = 0.5)  # Center the title
) +
guides(color = guide_legend(title = ""))
plot <- plot +
geom_vline(xintercept = as.Date("2005-12-01"), color = "black")
# Shading the area from 2006 onwards
# plot <- plot +
#   geom_rect(
#     xmin = as.Date("2006-01-01"),
#     xmax = max(plotforc3m$ACM.Date),
#     ymin = -Inf,
#     ymax = Inf,
#     fill = "gray",
#     alpha = 0.01,
#     color = NA  # Set the outline color to match the fill color
#   )
print(plot)
ggsave("figures/fig_fitforc_3m.pdf", plot = plot,width = 5.72, height = 3.50)
