Model Check - RSLN
Perform a check of the model to gain comfort. In this simplified example a call option is priced using the simulations as well as using the analytical Black-Scholes formula. The convergence will then be demonstrated.
Focus on call option with strike K. In this case, set K = S.
mu1 <- r #one-year swap rate
mu2 <- r #one-year swap rate
p12 <- 0 #probability of transitioning from state 1 to state 2
p21 <- 1 #probability of transitioning from state 2 to state 1
vol1 <- implied_vol
vol2 <- implied_vol
prices<- SimulateRSLNPrices(mu1, mu2, vol1, vol2, p12, p21, Nproj, dt, Nsim, calendar, S)
K <- S
payoffs <- as.vector(coredata(prices)[nrow(coredata(prices)),]) - K
payoffs <- ifelse(payoffs<0, 0, payoffs)
pv_payoffs <- payoffs*exp(-r)
mean(pv_payoffs)
## [1] 159.6744
bsp <- BS(S, S, r, 0, implied_vol, 1)
bsp
## [1] 160.1287
pv_payoff_converge <- cumsum(pv_payoffs)/index(pv_payoffs)
pv_payoff_converge.df <- data.frame(cbind(index(pv_payoffs),pv_payoff_converge))
pv_payoff_converge_pct.df <- data.frame(cbind(index(pv_payoffs),(pv_payoff_converge-bsp)/bsp))
colnames(pv_payoff_converge.df) <- c("Simulations", "AveragePresentValue")
colnames(pv_payoff_converge_pct.df) <- c("Simulations", "PctDifference")
pv_payoff_converge_pct.df <- subset(pv_payoff_converge_pct.df, Simulations>1000)
Model Check - Heston-Nandi GARCH(1,1)
Like Black-Scholes, the Heston-Nandi GARCH(1,1) model has analytical solutions for simple options. The call option examined above has a Heston-Nandi GARCH(1,1) price of 164.2191884. Below we see that the Heston-Nandi GARCH(1,1) simulations converge to the analytical solution and compare well with Black-Scholes.
payoffs <- as.vector(coredata(pricesHNG)[nrow(coredata(pricesHNG)),]) - K
payoffs <- ifelse(payoffs<0, 0, payoffs)
pv_payoffs <- payoffs*exp(-r)
mean(pv_payoffs)
## [1] 164.1752
hngo <- HNGOption(model=hngf$model, TypeFlag = "c", S=S, X=S, Time.inDays = Nproj, r.daily=log(1+r)/Nproj)
hngo
##
## Call:
## HNGOption(TypeFlag = "c", model = hngf$model, S = S, X = S, Time.inDays = Nproj,
## r.daily = log(1 + r)/Nproj)
##
## Option Price:
## 164.2192
pv_payoff_converge <- cumsum(pv_payoffs)/index(pv_payoffs)
pv_payoff_converge.df <- data.frame(cbind(index(pv_payoffs),pv_payoff_converge))
pv_payoff_converge_pct.df <- data.frame(cbind(index(pv_payoffs),(pv_payoff_converge-hngo$price)/hngo$price))
colnames(pv_payoff_converge.df) <- c("Simulations", "AveragePresentValue")
colnames(pv_payoff_converge_pct.df) <- c("Simulations", "PctDifference")
pv_payoff_converge_pct.df <- subset(pv_payoff_converge_pct.df, Simulations>1000)
Main Functions
SimulateRSLNPrices <- function(mu1, mu2, vol1, vol2, p12,p21, nproj, dt, nsim, calendar, startprice){
pi1<- p21/(p12+p21)
pi2 <- p12/(p12+p21)
#Generate random numbers to simlute the current state
stateRnd <- runif(nproj*nsim)
#Based on state, pick the corresponding return and vol
mu_sim<- matrix(ifelse(stateRnd< pi1,mu1,mu2), ncol=nproj, nrow=nsim)
vol_sim<- matrix(ifelse(stateRnd< pi1,vol1,vol2), ncol=nproj, nrow=nsim)
#Simulate the future prices
U <- matrix(exp((mu_sim-vol_sim^2/2)*dt+vol_sim*sqrt(dt)*rnorm(nproj*nsim)), ncol=nproj, nrow=nsim)
U <- t(apply(U, 1, cumprod))
U <- startprice*U
write.csv(U, "SimulateRSLNPrices.csv")
tU <- t(U)
tU.df <- data.frame(tU)
tU.xts <- xts(tU.df,order.by=calendar)
return(tU.xts)
}
hngarchSim.RN <- function (model = list(lambda = 4, omega = 4 * 2e-04, alpha = 0.3 *
2e-04, beta = 0.3, gamma = 0, rf = 0), n = 1000, innov = NULL,
n.start = 100, start.innov = NULL, rand.gen = rnorm, ...)
{
if (is.null(innov))
innov = rand.gen(n, ...)
if (is.null(start.innov))
start.innov = rand.gen(n.start, ...)
lambda = model$lambda
omega = model$omega
alpha = model$alpha
beta = model$beta
gamma = model$gamma
rfr = model$rf/n
x = h = Z = Z.star = c(start.innov, innov)
gamma.star <- gamma + lambda +.5
lambda = -.5
nt = n.start + n
h[1] = (omega + alpha)/(1 - alpha * gamma * gamma - beta)
x[1] = rfr + lambda * h[1] + sqrt(h[1]) * Z[1]
#Z.star[1] = Z[1] + (lambda+.5)*sqrt(h[1])
for (i in 2:nt) {
h[i] = omega + alpha * (Z[i - 1] - gamma.star * sqrt(h[i - 1]))^2 + beta * h[i - 1]
x[i] = rfr + lambda * h[i] + sqrt(h[i]) * Z[i]
}
x = x[-(1:n.start)]
x
print(mean(h))
print(mean(Z.star))
return(x)
}
SimulateHNGPrices <- function(model, nproj, nsim, calendar, startprice){
U <- create.matrix(nsim, nproj)
for(i in 1:nsim){
pricesHNG <- hngarchSim.RN(model=model, n=nproj)
U[i,]<- exp(pricesHNG)
}
U <- t(apply(U,1,cumprod))
U <- startprice*U
write.csv(U, "SimulateHNGPrices.csv")
tU <- t(U)
tU.df <- data.frame(tU)
tU.xts <- xts(tU.df,order.by=calendar)
return(tU.xts)
}
ApplyVolControl <- function(valdate, prices, riskfree, startvalue, nsims, ntrading, voldays1, voldays2, targetlevel, maxleverage, daylag, calendar, write.intermediate=FALSE){
valdate_loc <- match(valdate, index(prices))
start_price_hist_date <- index(prices)[valdate_loc-max(voldays1, voldays2)-daylag]
enddate_loc <- match(max(index(prices)), index(prices))
end_price_hist_date <- index(prices)[enddate_loc-daylag-1]
prices_short <- prices[paste0(as.character(start_price_hist_date),"::",as.character(end_price_hist_date))]
prices_shorter <- prices[paste0(as.character(valdate),"::")]
max_vol <- create.matrix(ntrading, nsims)
leverage <- create.matrix(ntrading, nsims)
max_leverage <- create.matrix(ntrading, nsims)
daily_index_returns <- create.matrix(ntrading, nsims)
daily_risk_control_returns <- create.matrix(ntrading, nsims)
final_levels <- create.matrix(ntrading, nsims)
risk_control_final_levels <- create.matrix(ntrading, nsims)
price_value <- create.matrix(1, nsims)
for(j in 1:nsims){
v_all_20 <- runSD(ROC(prices_short[,j]),n=voldays1, sample=FALSE)*sqrt(252)
v_all_40 <- runSD(ROC(prices_short[,j]),n=voldays2, sample=FALSE)*sqrt(252)
v_all <- na.omit(merge(v_all_20, v_all_40, all=FALSE))
check <- nrow(v_all)
max_vol[,j]<- apply(v_all, 1, max)
leverage[,j]<- target_level/max_vol[,j]
leverage[,j] <- matrix(sapply(leverage[,j], function(x) min(x, maxleverage)), nrow=ntrading, ncol=1)
pr <- periodReturn(prices_shorter[,j], period="daily", type="log", leading=FALSE)
tmp <- na.omit(as.matrix(coredata(pr)))
daily_index_returns[,j]<- tmp[,1]
for(i in 1:ntrading){
if(i ==1){
max_leverage[i,j] <- leverage[i,j]
}else{
if(max_leverage[i-1, j]-leverage[i,j] >max_change){
max_leverage[i,j] <- max_leverage[i-1, j] - max_change
}else{
if(leverage[i,j]-max_leverage[i-1, j]>max_change){
max_leverage[i,j] <- max_leverage[i-1,j] + max_change
}else{
max_leverage[i,j]<-leverage[i,j]
}
}
}
}
#print(j)
}
# Assumes that non-risky portion earns the swap rate
daily_risk_control_returns <- max_leverage*daily_index_returns + (1-max_leverage)*riskfree/ntrading
risk_control_final_levels <- exp(apply(daily_risk_control_returns, 2, cumsum))*startvalue
result <- mean(risk_control_final_levels[nrow(risk_control_final_levels),])*exp(-r)
if(write.intermediate){
print("Writing output")
write.zoo.transpose(prices, "prices.csv")
write.zoo.transpose(prices_short, "prices_short.csv")
write.csv(t(max_vol), "max_vol.csv")
write.csv(t(leverage), "leverage.csv")
write.csv(t(max_leverage), "max_leverage.csv")
write.csv(t(daily_risk_control_returns), "daily_risk_control_returns.csv")
write.csv(t(risk_control_final_levels), "risk_control_final_levels.csv")
write.csv(result, "pv.csv")
}
risk_control_final_levels <- xts(risk_control_final_levels, order.by = calendar)
return(risk_control_final_levels)
}
Helper Functions
getFred<-function(id, freq){
json_file <- paste0("https://api.stlouisfed.org/fred/series/observations?series_id=", id ,"&frequency=", freq , "&aggregation_method=eop&file_type=json&api_key=29fb098284eda4e97e4678353c8f25f2")
json_data <- fromJSON(txt=json_file)
myDat <- suppressWarnings(xts(as.numeric(json_data$observations$value), order.by=as.Date(json_data$observations$date)))
myDat <- na.omit(myDat)
colnames(myDat)<- id
assign(id, myDat, .GlobalEnv)
}
mult.xts <- function (x, y)
{
xmat <- as.matrix(x)
ymat <- as.matrix(y)
result <- c(rep(NA, NAs), result)
reclass(result, x)
}
create.matrix <- function(i, j) {
x <- matrix()
length(x) <- i*j
dim(x) <- c(i,j)
x
}
combine.past.proj <- function(past, proj, nsim){
tmp <- create.matrix(nrow(past), nsim)
for(i in 1:nsim){
tmp[,i] <- coredata(past)
}
tmp <- xts(tmp, order.by=index(past))
return(rbind(tmp, proj))
}
write.zoo.transpose<-function(x, file){
df<- data.frame(t(coredata(x)))
colnames(df)<-index(x)
write.csv(df,file=file)
}
