############# Review & Bootstrap: Estimation and C.I. (Code #2 - Used in Lecture 4) #############

############# Review: Hypothesis Testing  for IBM mean returns (capital gains) #############
SFX_da <- read.csv("https://www.bauer.uh.edu/rsusmel/4397/Stocks_FX_1973.csv", head=TRUE, sep=",")
names(SFX_da)
summary(SFX_da)

# Extract variables from imported data
x_ibm <- SFX_da$IBM			# extract IBM price data
x_RF <- SFX_da$RF			# extract Risk-free rate (in %)

# Define log returns & adjust size of variables accordingly
T <- length(x_ibm)			# sample size
lr_ibm <- log(x_ibm[-1]/x_ibm[-T])	# create IBM log returns (in decimal returns)

# T-test: H0: Mu_IBM = .10
alpha <- .10
freq <- 12 				# monthly frequency
mu_0 <- .05/freq
x_bar <- mean(lr_ibm)
t_test <- (x_bar - mu_0)/(sd(lr_ibm)/sqrt(T))
t_test

p_val <- (1 - pnorm(abs(t_test))) * 2		# pvalue for t_test (adjust b/c two sided test)
p_val


############# Review: (1-alpha)% C.I. for IBM mean returns #############
# 95% C.I.
z_alpha <- qnorm(1-alpha/2)				# z-value at alpha =.05 (2-sided test)
CI_lb <- x_bar - z_alpha * (sd(lr_ibm)/sqrt(T))		# Lower bound
CI_ub <- x_bar + z_alpha * (sd(lr_ibm)/sqrt(T))		# Upper bound
CI_lb
CI_ub

ci_95 <- c(CI_lb,  CI_ub)
ci_95

ci_95 * 12					# Annualized C.I.


############# Review: ERP Estimation and C.I. (with Shiller's historical dataset) #############

Sh_da <- read.csv("https://www.bauer.uh.edu/rsusmel/4397/Shiller_data.csv", head=TRUE, sep=",")
names(Sh_da)
SP <- Sh_da$P					# Extract P = S&P500 series
D <- Sh_da$D					# Extract D = Dividend S&P500 series
i <- Sh_da$Long_i				# Extract Long_i = 10-year interest rate series
T <- length(SP)
lr <- log((SP[-1]+D[-1]/12)/SP[-T])		# Define log returns
rf <- mean(i)/100/12
rf
T <- length (lr)


exc_ret <- lr - rf
var_exc_ret <- var(exc_ret)

erp <- mean(exc_ret)  				# mean monthly ERP
erp
erp * 12   				  	# annual ERP

# t-test
t_hat <- (erp - 0)/sqrt(var_exc_ret/T)
t_hat 

# 98% C.I.
z_98 <- 2.33					#z-value at alpha =.02
CI_lb <- erp - z_98 * sqrt(var_exc_ret/T)	# Lower bound
CI_ub <- erp + z_98 * sqrt(var_exc_ret/T)	# Upper bound
CI_lb
CI_ub

ci_98 <- c(erp - z_98 * sqrt(var_exc_ret/T),  erp + z_98 * sqrt(var_exc_ret/T))
ci_98

ci_98 * 12					# Annualized C.I.

## Note by looking at 98% C.I., 1% is outside range. This is Mehra & Prescott point, ERP is "too high"


############# Bootstrap: C.I. for return variance #############
sim_size <- 1000

## Bootstrap using boot ("perc") and defining function of statistic of interest
install.packages("boot")
library(boot)

# function to obtain the variance from the data
var_p <- function(data, i) {
			d <- data[i]
			return(var(d))
}

boot.samps <- boot(data=lr, statistic=var_p, R=sim_size)	# resampling and theta* estimation

## Percentile method using boot (Default is a 95% CI -i.e, conf = 0.95)
boot.ci(boot.samps, type = "perc")		# You can set conf (say, conf = 0.99)

## boot.ci(boot.samps, conf=.99, type = "perc")	

# Check CI bounds by sorting 1,000 variances (in boot.samps$t)
new <- sort(boot.samps$t)
sqrt(new[25]) 			 		# CI's Lower Bound 
sqrt(new[975]) 			 		# CI's Upper Bound 

hist(boot.samps$t, breaks=12, main="Histogram for Boostrapped Variances", xlab="Bootstrapped Variances")


## Empirical Bootstrap ("basic") using boot
boot.samps <- boot(data=lr, statistic=var_p, R=sim_size)	# resampling and theta* estimation
boot.ci(boot.samps, type = "basic")		# boot.ci computes the CI.

#Check bounds of Empirical Boostrap
lr_var <- var(lr)
q_star <- boot.samps$t - lr_var 		# q* = theta* −  theta^
q_ad <- sort(q_star)				# sort q*
lr_var - q_ad[975]				# CI's Lower Bound
lr_var - q_ad[25]			 	# CI's Upper Bound 			

sqrt(lr_var - q_ad[975])
sqrt(lr_var - q_ad[25])


#### Parametric Bootstrap assuming a Normal Dist
library(resample)
lr_var						# s2
lr_sd <- sqrt(lr_var)
m_lr <- mean(lr)				# mean
m_lr

x <- rnorm(T*sim_size, mean=m_lr, sd=lr_sd)		# generate normal data (mean & sd from data)
boot_par_sample <- matrix(x, nrow=T, ncol=sim_size)	# organize simulated data
boots_vars  <- colVars(boot_par_sample)			# compute variances
q_star <- boots_vars - lr_var				# center distribution around lr_var

q <- quantile(q_star, c(0.025, 0.975))
ci <- lr_var - c(q[2], q[1]) 
ci 
sqrt(ci)


# Sample SE of Var of bootstrap (used to gauge accuracy of bootstrap)
var(boots_vars)
sqrt(var(boots_vars))



#### Empirical C.I. for Sharpe Ratio
sr <- (mean(lr) - rf)/sd(lr)
sr				# Interpretation: a 1% in risk (SD) increases excess return by sr

### function to obtain the Sharpe Ratio from the data
sr_p <- function(data, i) {
			d <- data[i]
sr <- (mean(d) - rf)/sd(d)
			return(sr)
}

boot.samps <- boot(data=lr, statistic=sr_p, R=sim_size)	# resampling and θ* estimation
boot.ci(boot.samps, type = "basic")		# boot.ci computes the CI.

# Check Bounds
q_star <- boot.samps$t - sr			# q* = theta* -  theta^
q_ad <- sort(q_star)				# sort q*
sr- q_ad[975]					# CI's Lower Bound
sr - q_ad[25]			 		# CI's Upper Bound 			


# Percentile method
new <- sort(boot.samps$t)
new[25] 			 		# CI's Lower Bound 
new[975] 			 		# CI's Upper Bound 

hist(boot.samps$t, breaks=12, main="Histogram for Boostrapped Variances", xlab="Bootstrapped Variances")

#### VaR for FX Transaction

chfusd <- read.csv("http://www.bauer.uh.edu/rsusmel/4386/chfusd.csv",sep=",")  # Data
S <- chfusd$CHF_USD			# Extract CHF_USD column of the data
T <- length(S)				# Check total T (1971:1 to 2017:1)
Tstart <- 229				# Start of sample period: 1990:1
SP <- S[Tstart: T]			# FX Rate from relevante period (1990:1 on)
T <- length(SP)				# Number of observations 1990:1 to 2017:1
Val <- 1000000				# Value of transaction in FC (in M)
S_0 <- SP[T]				# FX Rate at T (Today's S_t)
e_f <- log(SP[-1]/SP[-T])		# Long changes in S_t
T_s <- length(e_f)			# Length of e_f (lost one observation)


## VaR assuming a Normal
TE <- Val *  S_0			# Transaction Exposure ("Amount Exposed")
mu_ef <- mean(e_f)			# Mean change in exchange rate (CHF/USD)
s_ef <- sd(e_f) 			# SD of change in exchange rate (CHF/USD)
mu_ef
s_ef
z_alpha2 <- qnorm(1-alpha/2)		# If alpha=05 => z_alpha2 = z_975 - 1.96
CI_lb <- mu_ef - z_alpha2 * s_ef	# Lower bound
CI_ub <- mu_ef + z_alpha2 * s_ef	# Upper bound
CI_lb
CI_ub

TE_lb <- TE * (1 + CI_lb)		# VaR (1-alpha/2) = Maximum loss within C.I.
TE_lb


## VaR using a Bootstrap
alpha = .05				# Specify alpha level for VaR(1-alpha/2)
T_s_low <- round(T_s * alpha/2)		# Obs corresponding to alpha/2*T_s
TE_o <- Val* S_0 * (1 + e_f)		# calculate TE for all e_f in sample ("Original TE" values 
STE_o <- sort(TE_o)			# sort Original TEs
VaR_o <- STE_o[T_s_low] 		# Pick obs in the 97.5th percentile ("Original VaR")

# function to obtain VaR from the data
varisk <- function(data, i) {
	d <-data[i]
  TE <- Val* S_0 * (1+d)		# calculate R TE values 
  STE <- sort(TE)			# sort TE
  VaR <- STE[T_s_low]			# Pick obs in the 97.5th percentile (975th obs if T_s = 1,000)
return(VaR)
}

library(boot)
sim_size <- 1000
conf_alpha <- 1 - alpha
boot.samps <- boot(data=e_f, statistic=varisk, R=sim_size)

boot.ci(boot.samps, conf = conf_alpha, type = "basic")

mean(boot.samps$t)			# Boostrapeed VaR(1-alpha/2)

sd(boot.samps$t)

# VaR-mean = Maximun expected loss within the one side (1-alpha/2) C.I.:
mean(boot.samps$t) - TE			# Maximum loss within C.I.