############# Review of VaR & Regression with lm function (Code #3A - Used in Lecture 5) #############

#### 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

# First, we build a (1-alpha)% C.I. for e_f
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 (" the worst case scenario" for a receivable within the C.I.)
CI_ub <- mu_ef + z_alpha2 * s_ef	# Upper bound ("the good/positive" bound for a receivable)
CI_lb
CI_ub

# Second, we use the worst case scenario to compute the VaR = Amount exposed * (1 + worst change in C.I for e_f)
TE_lb <- TE * (1 + CI_lb)		# Lower bound TE = VaR (1-alpha/2) = Maximum loss within C.I.

# Var-mean = Maximum Expected loss 
TE_lb - TE  				# Maximum Expected loss (within the C.I.)


## 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)

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

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



#### Regression using R with lm fuction

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_Mkt_RF <- SFX_da$Mkt_RF		# extract Market excess returns (in %)
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)
Mkt_RF <- x_Mkt_RF[-1]/100		# Adjust sample size to ( T-1) by removing 1st obs 
RF <- x_RF[-1]/100			# Adjust sample size and use decimal returns.

# Scatter plot of the two varialbes
plot(lr_ibm, Mkt_RF, main="Scatterplot: IBM & Market",
   xlab="IBM returns ", ylab="Market returns ", pch=19)

# Regression using lm
ibm_x <- lr_ibm - RF			# IBM excess returns
fit_ibm_capm <- lm(ibm_x ~ Mkt_RF)	# lm (=linear model) package in R
summary(fit_ibm_capm)			# print lm results

# Scatter plot of the two varialbes with regression line
plot(lr_ibm, Mkt_RF, main="Scatterplot: IBM & Market",
   xlab="IBM returns ", ylab="Market returns ", pch=19)
abline(fit_ibm_capm, col="blue")


# Regression using lm with R data frame, with different sample sizes
data_full <- data.frame(ibm_x, Mkt_RF)

fit_ibm_capm_full <- lm(ibm_x ~ Mkt_RF, data = data_full)	# lm reg with full data
summary(fit_ibm_capm_full)			
mean(Mkt_RF) * 12

T0 <- 1
T1 <- 200			# August 1989
data_short <- data_full[T0:T1, ]

fit_ibm_capm_short <- lm(ibm_x ~ Mkt_RF, data = data_short)	# lm reg with shorter data
summary(fit_ibm_capm_short)
mean(Mkt_RF[T0:T1]) * 12

T0 <- 201
T1 <- 400			# May 2006
data_short <- data_full[T0:T1, ]

fit_ibm_capm_short <- lm(ibm_x ~ Mkt_RF, data = data_short)	# lm reg with shorter data
summary(fit_ibm_capm_short)
mean(Mkt_RF[T0:T1]) * 12

T0 <- 401
T1 <- T				# End of Sample (Dec 2024)
data_short <- data_full[T0:T1, ]

fit_ibm_capm_short <- lm(ibm_x ~ Mkt_RF, data = data_short)	# lm reg with shorter data
summary(fit_ibm_capm_short)
mean(Mkt_RF[T0:T1]) * 12

## Note: Beta coeff are "time-varying" - The beta for the whole period is an average of betas over different periods
 (0.779049 * 223 + 1.076962 * 200 + 0.790583 * 200)/T
fit_ibm_capm_full$coefficients