Model_simulation(4).R

# Simulating stock returns and bond returns based on GDP regimes


#**********************************************************************
#                          Notes                                   ####
#**********************************************************************
# A discussion of ARIMA transfer function models:
# https://robjhyndman.com/hyndsight/arimax/


# Note
# Multivariate Normal
# https://www.r-bloggers.com/simulating-from-the-bivariate-normal-distribution-in-r/




## Issues:
  # Simulate ARIMA process with initial values and regressors

#**********************************************************************
#                           Packages                               ####
#**********************************************************************

library(markovchain) # Markov chain object
library(MASS)        # multivariate normal generator, must be loaded before tidyverse, otherwise 'select' will be masked


library(tidyverse)
library(broom)
library(readxl)
library(magrittr)
library(stringr)
library(forcats)
library(grid)
library(gridExtra)
library(scales)
library(knitr)

# packages for time series modeling
library(astsa)    # companion package
library(TSA)      # companion package;  arimax: flexible transfer function model
library(tseries)  #
library(forecast) # Arima
library(MSwM)
library(TTR)

#library(MSBVAR)


# packages for ts
library(zoo)
library(xts)


library(timetk)
library(tidyquant)

library(lubridate)
library(feather)

library(psych) # describe
library(xlsx)
library(qqplotr)

# check tidyquant, timetk, sweep (broom ), tibbletime
# Intro to zoo  cran.r-project.org/web/packages/zoo/vignettes/zoo-quickref.pdf
# sweep: http://www.business-science.io/code-tools/2017/07/09/sweep-0-1-0.html

#**********************************************************************
#                     Global settings                              ####
#**********************************************************************
dir_data_raw <- "data_raw/"
dir_data_out <- "data_out/"
dir_techReport <- "techReport_out"

# NBER recession periods, post-WWII
recessionPeriods <- 
  matrix(c(
  	1953+2/4, 1954+2/4,
  	1957+3/4, 1958+2/4,
  	1960+2/4, 1961+1/4,
  	1969+4/4, 1970+4/4,
  	1973+4/4, 1975+1/4,
  	1980+1/4, 1980+3/4,
  	1981+3/4, 1982+4/4,
  	1990+3/4, 1991+1/4,
  	2001+1/4, 2001+4/4,
  	2007+4/4, 2009+2/4
  ) , ncol = 2, byrow = T) %>% 
  as.data.frame() %>% 
	rename(peak =   V1,
				 trough = V2) %>% 
	mutate(peak = peak - 1/4,
				 trough = trough - 1/4)

get_geoReturn <- function(x) prod(1 + x)^(1/length(x)) - 1

RIG.theme <- function() {
	theme(
		panel.grid.major.x = element_blank(),
		panel.grid.minor.x = element_blank(),
		panel.grid.minor.y = element_blank(),
		panel.grid.major.y = element_line(size = 0.5, color = "gray80"),
		plot.title = element_text(hjust = 0.5),
		plot.subtitle = element_text(hjust = 0.5),
		plot.caption = element_text(hjust = 0, size = 9)
	)
}

RIG.themeLite <- function() {
	theme(
		plot.title = element_text(hjust = 0.5),
		plot.subtitle = element_text(hjust = 0.5),
		plot.caption = element_text(hjust = 0, size = 9)
	)
}



#**********************************************************************
#                     Loading Data                                 ####
#**********************************************************************
# Loading saved data 
load(paste0(dir_data_out, "dataAll.RData"))



#**********************************************************************
#               Data  ####
#**********************************************************************

# Stock return and volatility with different data frequencies

# Stock returns:
  # log total return (price appreciation + dividend)
  # Equity premium calculated with different measures of risk free rate
# Frequencies:
  # 1. Monthly
  # 2. Quarterly
  # 3. Yearly

get_logReturn <- function(x){
	if(any(x <= 0, na.rm = TRUE)) stop("Nagative value(s)")
	log(x/lag(x))
}

Vars <- c("year", "month", "yearMon", 
					"TBill3m_FRED", 
					"Tbond10y_FRED",
					"LCapStock_TRI",
					"LCapStock_CAI",
					"CBond_TRI",
					"LTGBond_TRI",
					"Inflation_Index",
					"GDP_FRED")


fn <- function(df, year_range, rolling_width, freq){
	df_stock_m <- 
		df %>% 
		select(one_of(Vars)) %>% 
		filter(year %in% year_range) %>% 
		mutate(return_tot = (1 + get_logReturn(LCapStock_TRI))^freq - 1,
					 dl_gdp     = (1 + get_logReturn(GDP_FRED))^freq - 1,
					 dl_cbond    = (1 + get_logReturn(CBond_TRI))^freq - 1,
					 dl_gbond    = (1 + get_logReturn(LTGBond_TRI))^freq - 1,
					 
					 return_tot_o = get_logReturn(LCapStock_TRI),
					 dl_gdp_o     = get_logReturn(GDP_FRED),
					 dl_cbond_o   = get_logReturn(CBond_TRI),
					 dl_gbond_o   = get_logReturn(LTGBond_TRI),
					 
					 TBill3m_FRED =  TBill3m_FRED/100,
					 Tbond10y_FRED = Tbond10y_FRED/100,
					 ERP_3m     = return_tot - TBill3m_FRED,
					 ERP_10y    = return_tot - Tbond10y_FRED) %>% 
		mutate(sd_return    = rollapply(return_tot, rolling_width, sd,   align = "right", fill = NA),
					 mean_return  = rollapply(return_tot, rolling_width, mean, align = "right", fill = NA),
					 sd_ERP_3m    = rollapply(ERP_3m, rolling_width, sd,    align = "right", fill = NA),
					 mean_ERP_3m  = rollapply(ERP_3m, rolling_width, mean,  align = "right", fill = NA),
					 sd_ERP_10y   = rollapply(ERP_10y, rolling_width, sd,   align = "right", fill = NA),
					 mean_ERP_10y = rollapply(ERP_10y, rolling_width, mean, align = "right", fill = NA),
					 sd_cbond     = rollapply(dl_cbond, rolling_width, sd,   align = "right", fill = NA),
					 mean_cbond   = rollapply(dl_cbond, rolling_width, mean, align = "right", fill = NA),
					 sd_gbond     = rollapply(dl_gbond, rolling_width, sd,   align = "right", fill = NA),
					 mean_gbond   = rollapply(dl_gbond, rolling_width, mean, align = "right", fill = NA),
					 mean_gdp     = rollapply(dl_gdp, rolling_width, mean,  align = "right", fill = NA)
		)
}

df_stock_m <- fn(df_dataAll,   1953:2015, 12, 12)
df_stock_q <- fn(df_dataAll_q, 1953:2015, 12, 4)
df_stock_y <- fn(df_dataAll_y, 1953:2015, 5,  1)


# save data in feather format for python use
write_feather(df_stock_q, "data_out/df_stock_q.feather" )
df_stock_q$dl_gdp_o



#***********************************************************************************
#                    0. Overview the simulation procedure          ####
#***********************************************************************************

#' Overview of the simplified modeling approach:
#' 
#' 1. GDP (quarterly, core variable):
#'    - Markov-switching model (AR or RW-drift)
#'    - Transition matrix that determines regime-switching behavior
#'    - Regime-dependent mean/std/parameters
#'         - mean/std based on NBER definition
#'         - mean/std based on stock MS model
#'         - mean/std based on GDP   MS model
#'    - Simulated regimes and paths of GDP growth
#' 2. Total stock returns (for now, we may want to model ERP later in order to incorporate inflation)
#'    - Markov-switching model: RW-drift
#'    - Return regimes are highly aligned with GDP regimes
#'    - For simulation: different mean return and std in GDP different GDP regimes
#' 3. Total cbond returns 
#'    - Show regime-switching behavior, but very different from those for GDP and stock returns.
#'       - high return and high volatility regime 
#'       - low return and low volatility regime
#'       - May not be pure random walk within regimes, may need to model MSAR process. 
#'    - Correlation with stock changes wildly over time, with no obivous pattern. 
#'    - Simulation stragegy:
#'       - 1. modeled and simulated as a seperate MS-AR/DR process, calibrate 
#'            correlation with stock to historical value/assumed value by setting
#'            the correlation between error terms of stock model and bond model.
#'       - 2. 
#' 4. Return of Portfolio
#'    - Construct portfolio returns   
#


#' Summary of quaterly stock mean return and std 
#' 
#' 1. Based on NBER recession periods 1953 - 2015
#'            mean       std
#' expansion  0.032     0.0686
#' recession -0.014     0.119  
#' 
#' 
#' 2. Based on MS model of stock return 1954-2015
#'            mean       std
#' expansion:  0.0419  0.0547
#' recession: -0.0192  0.1126
#' 
#' 
#' 3. Based on MS RW-drift model of GDP (regimes imported from python program)
#' 
#' - estimate seed 123 (?a local optimum, but consistent with hamilton1989)
#'              mean       std
#' expansion:  0.032   0.0685
#' recession:  0.0049  0.1044
#' 
#' - estimate seed 127 (?global optimum, matches historical length of recessions best)
#'             mean       std
#' expansion: 0.028    0.0723
#' recession: 0.0027   0.122
#' 



#' Comparing NBER and model recession/expansion period length  
#' 
#' US recessions NBER:
#' 	recession       expansion
#' Starting 1953:      3.7Q            20.2Q
#' Starting 1980       3.7Q            23.7Q
#' 
#' 
#' Model recessions
#'                   recession       expansion
#' MSAR4               1.6             19.9 
#' MSAR2               1.7             20.9      
#' MSAR1               1.8             20.4
#' RW-drift            4.4             11.4
#' RW-drift            3.1             20.7 
#' Hamilton1989        4.1             16.2


#' Summary  
#' 
#' - Simulated quarterly returns
#'    - Simulated mean and annual returns are reasonable and generally consistent with historical data
#'    - Skewness and kurtosis are still much smaller than historical  
#'    - Compare to normal: 
#'    
#'    
#' - Simulated annual returns
#'    - Simulated mean annual return and standard deviation looks reasonable
#'    - Kurtosis of historical kurtosis is surprisingly low 
#'    - compare to normal"
#'    
#' - Note: 
#'    - qarterly returns non-normal and fat tail
#'    - It seems that annual returns are normal?


#'Type of returns
#'  Frequency:
#'   1. quarterly
#'   2. annual
#'  How returns are generated:
#'   1. historical data
#'   2. simulations using gdp regimes
#'     2.1 mean/std from historical regimes
#'     2.2 mean/std from gdp model regimes
#'     2.3 mean/std from stock model regimes
#'     2.4 other scenarios
#'   3.Normal distribution
#'   
#'Measures and plots
#' mean
#' std
#' skewness
#' kurtosis
#' 1%, 5% percentiles
#' 
#' histogram
#' kernal smoothed curve
#' qqplot
#'   




#***********************************************************************************
#                 1. Global parameters for simulation                    ####
#***********************************************************************************

nyear <- 30
nsim  <- 2000


#***********************************************************************************
#             2. Defining the 2-state Markov chain for GDP             ####
#***********************************************************************************

# GDP regimes:
gdpStates <- c('0',  # 0: Expansion
							 '1')  # 1: recession


# Transition probabilities of GDP 
p01 <- 0.048 # Expansion to recession
p10 <- 0.32  # Recession to expansion 
p00 <- 1 - p01
p11 <- 1 - p10

gdpMatrix <- matrix(c(p00, p01,
											p10, p11), 
										 byrow = TRUE, 
										 nrow = 2) 
gdpMatrix


# Defining markovchain object
mc_gdp <- new("markovchain", 
							states = gdpStates,
							transitionMatrix = gdpMatrix,
							name = 'gdp'
							)

mc_gdp %>% plot
summary(mc_gdp)

# Try simulation with markov chain
rmarkovchain(n = 4*30, object = mc_gdp, t0 = '0')
replicate(10, rmarkovchain(n = 4*30, object = mc_gdp, t0 = '1')) 




#***********************************************************************************
#             3. Paramters of gdp growth, stock and bond returns         ####
#***********************************************************************************

# # Stock: expansion
# s.mean_0 <- 0.032
# s.std_0  <- 0.069 
# 
# # Stock: recession
# s.mean_1 <- -0.014
# s.std_1  <- 0.119
# 
# # Bond
# b.mean  <- 0.016 #0.035/4
# b.std   <- 0.051 #0.04/2
# 
# # Correlation between error terms(deviation from mean) of stock and bond returns
# # (Calibration needed)
# # Goal is to find a rho that can make the correlation between the simulated
# # annual stock returns and annual bond returns (both calculated from quarterly returns)
# # match the target value
# rho <- 0.15 
# 
# # GDP expansion 
# gdp.mean_0 <- 0.0095
# gdp.std_0  <- 5.463e-05^0.5 # 0.007391211
# 
# # GDP recession
# gdp.mean_1 <- -0.0055
# gdp.std_1  <- 5.463e-05^0.5 # 0.007391211



simInputs_historical <- 
	list(
		nyear = 30,
		nsim  = 2000,
		
		# Markov Chain object for GDP
		mc_gdp = mc_gdp,
		
		# Stock: expansion
		s.mean_0 = 0.032,
		s.std_0  = 0.069 ,
		
		# Stock: recession
		s.mean_1 = -0.014,
		s.std_1  = 0.119,
		
		# Bond
		b.mean  = 0.016, #0.035/4
		b.std   = 0.051, #0.04/2
		
		# Correlation between error terms(deviation from mean) of stock and bond returns
		# (Calibration needed)
		# Goal is to find a rho that can make the correlation between the simulated
		# annual stock returns and annual bond returns (both calculated from quarterly returns)
		# match the target value
		rho = 0.15, 
		
		# GDP expansion 
		gdp.mean_0 = 0.0095,
		gdp.std_0  = 5.463e-05^0.5, # 0.007391211
		
		# GDP recession
		gdp.mean_1 = -0.0055,
		gdp.std_1  = 5.463e-05^0.5 # 0.007391211
		
	)
	

simInputs_forward <- 
	list(
		nyear = 30,
		nsim  = 2000,
		
		# Markov Chain object for GDP
		mc_gdp = mc_gdp,
		
		# Stock: expansion
		s.mean_0 = 0.032 - 0.00656,  #  0.02544
		s.std_0  = 0.069 ,
		
		# Stock: recession
		s.mean_1 = -0.014 - 0.00656, # -0.0205
		s.std_1  = 0.119,
		
		# Bond
		b.mean  = 0.009086, # 
		b.std   = 0.02, #  0.04/2
		
		# Correlation between error terms(deviation from mean) of stock and bond returns
		# (Calibration needed)
		# Goal is to find a rho that can make the correlation between the simulated
		# annual stock returns and annual bond returns (both calculated from quarterly returns)
		# match the target value
		rho = 0.15, 
		
		# GDP expansion 
		gdp.mean_0 = 0.0095 - 0.002786, # 0.006714
		gdp.std_0  = 5.463e-05^0.5,     # 0.007391211
		
		# GDP recession
		gdp.mean_1 = -0.0055 - 0.002786, # -0.008286
		gdp.std_1  = 5.463e-05^0.5       # 0.007391211
	)

# GDP 
  # Target: 1.019242^0.25 - 1 = 0.004776
  # simulated based on historical data:0.007562
  # adjustment: 0.007562 - 0.004776 = 0.002786

# Stock:
  # Target:  1.08145^0.25 - 1 =  0.01976854
  # current: 1.1095593^0.25 - 1 = 0.02633
  # adjustment: 0.00656146


#***********************************************************************************
#     5. Simulating quarterly GDP growth and asset returns with regime-switching   ####
#***********************************************************************************

# simInputs <- simInputs_historical
# sim_name  <- "sim_results_historical"

simInputs <- simInputs_forward
sim_name  <- "sim_results_forward"




set.seed(11)
{
	
# Generating GDP regimes	
sim_gdp_regimes <- replicate(simInputs$nsim, rmarkovchain(n = simInputs$nyear * 4, object = simInputs$mc_gdp, t0 = '0') %>% as.numeric())

	
# Generating of stock and bond returns (approach 1)
{
 # Mean vector and covariance matrix 
 # (Assuming the correlation between bond and stock returns does not change over regimes)
# mu_0 <- c(s.mean_0, b.mean)
# mu_1 <- c(s.mean_1, b.mean)
# 
# sigma_0 <- matrix(c(s.std_0^2,            s.std_0*b.std*rho,
# 							 		  s.std_0*b.std*rho,    b.std^2), 
# 								    2)
# sigma_1 <- matrix(c(s.std_1^2,            s.std_1*b.std*rho,
# 										s.std_1*b.std*rho,    b.std^2), 
# 									  2)
# mu;sigma_0;sigma_1
# 
# sim_errorTerms_0 <- replicate(nsim, mvrnorm(nyear, mu = mu_0, Sigma = sigma_0)) 
# sim_errorTerms_1 <- replicate(nsim, mvrnorm(nyear, mu = mu_1, Sigma = sigma_1)) 
# 
# sim_stockreturn_0 <- sim_errorTerms_0[,1,]
# sim_stockreturn_1 <- sim_errorTerms_1[,1,]
# 
# sim_bondreturn_0 <- sim_errorTerms_0[,2,]
# sim_bondreturn_1 <- sim_errorTerms_1[,2,]
}
	
# Generating of stock and bond returns (approach 2, based on standard multivariate normal distribution)

sim_errorTerms_stdNormal <- replicate(simInputs$nsim, mvrnorm(simInputs$nyear*4, mu = c(0, 0), Sigma = matrix(c(1, simInputs$rho, simInputs$rho, 1), 2))) 

sim_stockreturn_0 <- sim_errorTerms_stdNormal[,1,]*simInputs$s.std_0 + simInputs$s.mean_0
sim_stockreturn_1 <- sim_errorTerms_stdNormal[,1,]*simInputs$s.std_1 + simInputs$s.mean_1

sim_bondreturn <- sim_errorTerms_stdNormal[,2,]*simInputs$b.std + simInputs$b.mean

# corr.test(cbind(as.vector(sim_bondreturn), as.vector(sim_stockreturn_1))) %>% print(short = F)
# sim_stockreturn_0 <- replicate(nsim, rnorm(nyear*4, s.mean_0, s.std_0))
# sim_stockreturn_1 <- replicate(nsim, rnorm(nyear*4, s.mean_1, s.std_1))

#  Simulating GDP growth
sim_gdp_0 <- replicate(simInputs$nsim, rnorm(simInputs$nyear*4, simInputs$gdp.mean_0, simInputs$gdp.std_0))
sim_gdp_1 <- replicate(simInputs$nsim, rnorm(simInputs$nyear*4, simInputs$gdp.mean_1, simInputs$gdp.std_1))

# Compute stock returns with regime-switching
sim_stockreturn <- (sim_gdp_regimes == 0) * sim_stockreturn_0 + (sim_gdp_regimes == 1) * sim_stockreturn_1
sim_gdp_growth  <- (sim_gdp_regimes == 0) * sim_gdp_0 +         (sim_gdp_regimes == 1) * sim_gdp_1
}


# check correlation between regime-switching stock returns and single-regime bond returns
cor(cbind(as.vector(sim_bondreturn), as.vector(sim_stockreturn))) 


# sim_stockreturn
# sim_bondreturn
# sim_gdp_growth
# 
# sim_gdp_regimes
# 
# sim_stockreturn_0 
# sim_stockreturn_1 
# 
# sim_gdp_0 %>% mean 
# sim_gdp_1 %>% mean



#' TO what extent the simulation approach is capable of generating economic and investment return
#' scenarios with statistical characteristics that are similar to historical data.

#sim_gdp_regimes %>% dim()


#***********************************************************************************
#           6. Converting results to annual data and save                         ####
#***********************************************************************************

## Convert quarterly returns / growth into annual values (quarterly returns compounded within a year)

df_sim_stockreturn_q <- 
	as.data.frame(sim_stockreturn) %>% 
	mutate(year = rep(1:nyear, each = 4))
df_sim_stockreturn_y <-  
	df_sim_stockreturn_q %>% 
	gather(sim, return, -year) %>% 
	group_by(sim, year) %>% 
	summarise(return_y = prod(1+return) - 1 )


df_sim_bondreturn_q <- 
	as.data.frame(sim_bondreturn) %>% 
	mutate(year = rep(1:nyear, each = 4))
df_sim_bondreturn_y <-  
	df_sim_bondreturn_q %>% 
	gather(sim, return, -year) %>% 
	group_by(sim, year) %>% 
	summarise(return_y = prod(1+return) - 1 )

df_sim_gdp_q <- 
	as.data.frame(sim_gdp_growth) %>% 
	mutate(year = rep(1:nyear, each = 4))
df_sim_gdp_y <-  
	df_sim_gdp_q %>% 
	gather(sim, return, -year) %>% 
	group_by(sim, year) %>% 
	summarise(return_y = prod(1+return) - 1 )

# regimes
df_sim_gdp_regimes_q <- 
	as.data.frame(sim_gdp_regimes) %>% 
	mutate(year = rep(1:nyear, each = 4))


df_sim_gdp_regimes_y <- 
	as.data.frame(df_sim_gdp_regimes_q) %>% 
	gather(sim, regime, -year) %>% 
	group_by(sim, year) %>% 
	summarise(recession_nqtr = sum(regime))


sim_results <- list(
	df_sim_stockreturn_q = df_sim_stockreturn_q,
	df_sim_bondreturn_q  = df_sim_bondreturn_q,
	df_sim_gdp_q         = df_sim_gdp_q,
	df_sim_gdp_regimes_q = df_sim_gdp_regimes_q,
	
	
	df_sim_stockreturn_y= df_sim_stockreturn_y,
	df_sim_bondreturn_y = df_sim_bondreturn_y,
	df_sim_gdp_y        = df_sim_gdp_y,
	df_sim_gdp_regimes_y= df_sim_gdp_regimes_y
)


assign(sim_name, sim_results)


# save(df_sim_stockreturn_q,
# 	   df_sim_bondreturn_q,
# 	   df_sim_gdp_q,
# 	   df_sim_gdp_regimes_q,
# 	   
# 	   
# 	   df_sim_stockreturn_y,
# 	   df_sim_bondreturn_y,
# 	   df_sim_gdp_y,
# 	   df_sim_gdp_regimes_y,
save(sim_results,
		 file = paste0("policyBrief_out/", sim_name, ".RData"))

sim_name


#***********************************************************************************
#          7. Examine recession regime: number and length in 30-year period        ####
#***********************************************************************************

# calculate length of regimes in each simulation
get_regimeLength <- function(x, rgm = 1){
	
	#x <- sim_gdp_regimes[,2]
	y <- numeric(length(x))
	y[1] <- 1
	
	for(i in 2:length(x)){
		if(x[i] != x[i-1]) y[i] <- y[i-1] + 1 else y[i] <- y[i-1] 
	}
	
	df <- data.frame(idx = y, regime = x)
	df %<>% filter(regime == rgm) %>% 
		group_by(idx) %>% 
		summarise(regime_len = n())
	df
	
	r <- df$regime_len
	r
}

summary_gdpRegime_0 <- sapply(as.data.frame(sim_gdp_regimes), get_regimeLength, rgm = 0)
summary_gdpRegime_1 <- sapply(as.data.frame(sim_gdp_regimes), get_regimeLength, rgm = 1)

# Number of recessions 
sapply(summary_gdpRegime_1, length) %>% mean  # ~5 recessions in 30 years (many are very short)

# Expected length of regimes
sapply(summary_gdpRegime_1, mean) %>% mean(., na.rm = TRUE)  # average length is ~3.1 quarters. 
sapply(summary_gdpRegime_0, mean) %>% mean(., na.rm = TRUE)  # average length is ~20  quarters. 
# all as expected

# distribution of recession length (all sims together)
x <- data.frame(n.rec = unlist(summary_gdpRegime_1))
x %>% count(n.rec)
nrow(x)	
unlist(summary_gdpRegime_1) %>% hist()
# about 1/3 "recessions" only last for 1 quarter

# distribution of average recession length across sims
sapply(summary_gdpRegime_1, mean) %>% hist()

# distribution of number of recessions across sims
sapply(summary_gdpRegime_1, length) %>% hist()




#***********************************************************************************
#           8. Examine simulation results                        ####
#***********************************************************************************

# ## Check annual returns and growth
# 
# annual stock returns, all sims
df_sim_stockreturn_y$return_y	%>% mean # 10.8% mean annual return (quarterly return compounded)
df_sim_stockreturn_y$return_y	%>% sd   # 17.2% std
df_sim_stockreturn_y$return_y	%>% describe  # ~0.57 kurtosis, higher than historical returns
df_sim_stockreturn_y$return_y %>% hist(seq(-0.65, 1.17, 0.02))

df_sim_stockreturn_y %>% 
	group_by(sim) %>% 
	summarise(geoReturn = get_geoReturn(return_y)) %>% 
	summarise(geoReturn_mean = mean(geoReturn)) # 6.8%
df_sim_stockreturn_y$return_y	%>% sd   # 17.2% std

# annual bond returns, all sims
df_sim_bondreturn_y$return_y	%>% mean # ~3.5% mean annual return (quarterly return compounded)
df_sim_bondreturn_y$return_y	%>% sd   # ~4.1% std
df_sim_bondreturn_y$return_y	%>% describe  # very small skewness and kurtosis
df_sim_bondreturn_y$return_y %>%  hist(seq(-0.25, 0.25, 0.01))


df_sim_bondreturn_y %>% 
	group_by(sim) %>% 
	summarise(geoReturn = get_geoReturn(return_y)) %>% 
	summarise(geoReturn_mean = mean(geoReturn)) # 3.6%
df_sim_bondreturn_y$return_y	%>% sd   # 4.1% std



# annual gdp returns, all sims
df_sim_gdp_y$return_y	%>% mean # ~3.05% mean annual return (quarterly return compounded)
df_sim_gdp_y$return_y	%>% sd   # ~2.2% std (larger than std of shocks due to regime-switching)
df_sim_gdp_y$return_y	%>% describe  # heavily skewed, 0.8 kurtosis
df_sim_gdp_y$return_y %>% hist
# 

df_sim_gdp_y %>% 
	group_by(sim) %>% 
	summarise(geoReturn = get_geoReturn(return_y)) %>% 
	summarise(geoReturn_mean = mean(geoReturn)) # 1.9%
df_sim_gdp_y$return_y	%>% sd   # 2.2% std


cor(df_sim_stockreturn_y$return_y,df_sim_bondreturn_y$return_y)


df_sim_gdp_y





## Compared with historical distribution

## qqplot for GDP growth
library(qqplotr)

nsim_plot <- 50

df_qqplot_sim <-
bind_rows(
	df_stock_y %>% select(value = dl_gdp_o)    %>% mutate(type = "Historical", var = 'GDP growth' ),
	(df_sim_gdp_y %>% select(value = return_y) %>% mutate(type = 'Simulated',  var = 'GDP growth'))[1:(nsim_plot*nyear),],  # plot based on 50 simulations

	df_stock_y %>% select(value = return_tot_o ) %>% mutate(type = "Historical", var = 'Stock return'),
	(df_sim_stockreturn_y %>% select(value = return_y) %>% mutate(type = 'Simulated',   var = 'Stock return'))[1:(nsim_plot*nyear),],  # plot based on 50 simulations

	df_stock_y           %>% select(value = dl_gbond_o) %>% mutate(type = "Historical", var = 'Bond return'),
	(df_sim_bondreturn_y %>% select(value = return_y)   %>% mutate(type = 'Simulated',  var = 'Bond return'))[1:(nsim_plot*nyear),]  # plot based on 50 simulations
	)
df_qqplot_sim %>% head


fig_qqplot_simGDP <-
	df_qqplot_sim %>% filter(var == 'GDP growth') %>%
	ggplot(aes(sample = value)) + facet_wrap(var~type, scales = 'fixed') + theme_bw() + RIG.themeLite()+
	stat_qq_point(size = 1) +
	stat_qq_line() +
	#stat_qq_band(alpha = 0.5, con = 0.95, bandType = "boot") +
	labs(x = "Theoretical Quantiles", y = "Sample Quantiles",
			 title = 'Comparing Q-Q plots of historical annual GDP growth and simulated annual GDP growth')
fig_qqplot_simGDP


fig_qqplot_simStock <-
	df_qqplot_sim %>% filter(var == 'Stock return') %>%
	ggplot(aes(sample = value)) + facet_wrap(var~type, scales = 'fixed') + theme_bw() + RIG.themeLite()+
	stat_qq_point(size = 1) +
	stat_qq_line() +
	#stat_qq_band(alpha = 0.5, con = 0.95, bandType = "boot") +
	labs(x = "Theoretical Quantiles", y = "Sample Quantiles",
			 title = 'Comparing Q-Q plots of historical annual stock return and simulated annual stock return')
fig_qqplot_simStock


fig_qqplot_simBond <-
	df_qqplot_sim %>% filter(var == 'Bond return') %>%
	ggplot(aes(sample = value)) + facet_wrap(var~type, scales = 'fixed') + theme_bw() + RIG.themeLite()+
	stat_qq_point(size = 1) +
	stat_qq_line() +
	#stat_qq_band(alpha = 0.5, con = 0.95, bandType = "boot") +
	labs(x = "Theoretical Quantiles", y = "Sample Quantiles",
			 title = 'Comparing Q-Q plots of historical annual bond return and simulated annual bond return')
fig_qqplot_simBond



ggsave(paste0(dir_techReport, "fig_qqplot_simGDP.png"),   fig_qqplot_simGDP,   width = 10*0.95, height = 5.5*0.95)
ggsave(paste0(dir_techReport, "fig_qqplot_simStock.png"), fig_qqplot_simStock, width = 10*0.95, height = 5.5*0.95)
ggsave(paste0(dir_techReport, "fig_qqplot_simBond.png"),  fig_qqplot_simBond,  width = 10*0.95, height = 5.5*0.95)


df_sim_descriptive <-
df_qqplot_sim %>%
	filter(!is.na(value)) %>%
	group_by(var, type) %>%
	select(-sim) %>%
	do(describe(.$value))

df_sim_descriptive

write_csv(df_sim_descriptive, paste0(dir_techReport, 'descriptive_stats_sim.csv'))





(df_sim_gdp_y %>% select(GDP_growth = return_y) %>% mutate(type = 'Simulated'))[1:900,]$GDP_growth %>% describe



## Quantiles of GDP simulation

df_sim_gdp_y %>%
	group_by(sim) %>%
	summarise(avg = mean(return_y),
						std  = sd(return_y)) %>%
	summarise(avg_q25 = quantile(avg, 0.25),
						avg_q50 = quantile(avg, 0.50),
						avg_q75 = quantile(avg, 0.75),

						std_q25 = quantile(std, 0.25),
						std_q50 = quantile(std, 0.50),
						std_q75 = quantile(std, 0.75)
						)

# number of recessions

# Number of recessions
sapply(summary_gdpRegime_1, length) %>% quantile(0.25)   # ~5 recessions in 30 years (many are very short)
sapply(summary_gdpRegime_1, length) %>% quantile(0.50)   # ~5 recessions in 30 years (many are very short)
sapply(summary_gdpRegime_1, length) %>% quantile(0.75)   # ~5 recessions in 30 years (many are very short)

# Expected length of regimes
sapply(summary_gdpRegime_1, mean) %>% quantile(0.25, na.rm = TRUE)   # average length is ~3.1 quarters.
sapply(summary_gdpRegime_1, mean) %>% quantile(0.50, na.rm = TRUE)   # average length is ~3.1 quarters.
sapply(summary_gdpRegime_1, mean) %>% quantile(0.75, na.rm = TRUE)   # average length is ~3.1 quarters.

sapply(summary_gdpRegime_0, mean) %>% quantile(0.25, na.rm = TRUE)   # average length is ~3.1 quarters.
sapply(summary_gdpRegime_0, mean) %>% quantile(0.50, na.rm = TRUE)   # average length is ~3.1 quarters.
sapply(summary_gdpRegime_0, mean) %>% quantile(0.75, na.rm = TRUE)   # average length is ~3.1 quarters.



sapply(summary_gdpRegime_1, mean) %>% head

summary_gdpRegime_1 %>% head

sapply(summary_gdpRegime_0, mean) %>% mean(., na.rm = TRUE)  # average length is ~20  quarters.
# all as expected




## Stock
# qq-plot
qqnorm(rnorm(10000))
qqline(rnorm(10000))

qqnorm(df_sim_stockreturn_y$return_y[1:10000])
qqline(df_sim_stockreturn_y$return_y[1:10000])

df_stock_y$return_tot_o[-1] %>% qqnorm
df_stock_y$return_tot_o[-1] %>% qqline

# descriptive statistics
rnorm(10000) %>% describe
df_sim_stockreturn_y$return_y[1:10000] %>% describe
df_stock_y$return_tot_o[-1] %>% describe


df_stock_q$return_tot_o[-1] %>% qqnorm
df_stock_q$return_tot_o[-1] %>% qqline

# probability of -30% return or worse
pnorm(-0.3, 0.108, 0.172) # 0.8%
pnorm(-0.3, 0.067, 0.164)


(df_sim_stockreturn_y$return_y[1:10000]<= -0.3) %>% sum


df_sim_stockreturn_y$return_y[1:10000]


## bond

df_sim_bondreturn_y$return_y[1:10000] %>% qqnorm
df_sim_bondreturn_y$return_y[1:10000] %>% qqline

df_stock_y$dl_cbond[-1] %>% qqnorm
df_stock_y$dl_cbond[-1] %>% qqline


## gdp growth

df_sim_gdp_y$return_y[1:10000] %>% qqnorm
df_sim_gdp_y$return_y[1:10000] %>% qqline

df_stock_y$dl_gdp[-1] %>% qqnorm
df_stock_y$dl_gdp[-1] %>% qqline



df_stock_y





# simulated quarterly stock returns

sim_stockreturn %>% mean
sim_stockreturn %>% sd
sim_stockreturn %>% as.vector() %>% describe


sim_gdp_growth %>% mean
sim_gdp_growth %>% sd





df_sim_stockreturn_y$return_y	     %>% hist(27)
sim_stockreturn	%>% as.vector()    %>% hist(27)
rnorm(20000,0.02607146, 0.0790173) %>% hist(27)


rnorm(20000,0.02607146, 0.0790173) %>% quantile(0.10)
sim_stockreturn %>% quantile(0.10)


# historical data (quarterly)

#stock quarterly
df_stock_q$return_tot_o[-1] %>% mean #2.55%
df_stock_q$return_tot_o[-1] %>% sd   #7.9%
df_stock_q$return_tot_o[-1] %>% describe # skew -0.92, kurtosis 1.72
rnorm(2000,0.0255,0.0792)   %>% describe

# stock annual
df_stock_y$return_tot_o[-1] %>% mean  # 10.3%
df_stock_y$return_tot_o[-1] %>% sd    # 15.7%
df_stock_y$return_tot_o[-1] %>% describe() # skew -0.38, kurtosis 0.03
rnorm(20000,0.0255,0.0792)   %>% describe   # skew 0.03,  kur

df_stock_y$return_tot_o[-1] %>% plot(type = 'l')

df_stock_y$return_tot_o[-1] %>% hist(27)
df_stock_q$return_tot_o[-1] %>% hist(27)

df_stock_q$return_tot_o[-1] %>% qqnorm
df_stock_q$return_tot_o[-1] %>% qqline

df_stock_y$return_tot_o[-1] %>% qqnorm
df_stock_y$return_tot_o[-1] %>% qqline




df_stock_q$return_tot_o[-1] %>% quantile(0.10)



#cbond
df_stock_q$dl_cbond_o[-1] %>% mean
df_stock_q$dl_cbond_o[-1] %>% sd

#gbond
df_stock_q$dl_gbond_o[-1] %>% mean
df_stock_q$dl_gbond_o[-1] %>% sd

shapiro.test(df_stock_q$return_tot_o)





#***********************************************************************************
#                             Saving results                                ####
#***********************************************************************************

# save(
# df_sim_stockreturn_q,
# df_sim_bondreturn_q,
# df_sim_gdp_q,
# df_sim_gdp_regimes_q,
# 
# 
# df_sim_stockreturn_y,
# df_sim_bondreturn_y,
# df_sim_gdp_y,
# df_sim_gdp_regimes_y,
# file = "policyBrief/simulation_MS1.RData")
# 
# 
# load("techReport_out/simulation_MS1.RData")