R261_Homework 用歷史模擬法計算VaR值
AllenTsai
2015年11月9日
歷史模擬法：衍生性商品的應用
選取一段時間內即期匯率之價格
計算此段時間即期匯率之報酬, r1, r2,..r100
計算所有隔日即期匯率的可能 \({x = y (e^{r_i} - 1)}\)
計算本日衍生性商品的 Delte、Gamma
計算隔日衍生性商品之價格變動，並由小到大排序
設定信賴區間(假設 95%)
VaR = 累積百分比達到5%的價格
Black Sholes Formula
\({d_1}\)
\({\frac{log(S_0 / K) + (r -q + 0.5 * \sigma ^2 ) * t } {\sigma * \sqrt{t}} }\)
getD1 <- function(spot,k,rd,rf,sigma,t){
(log(spot/k)+(rd-rf+0.5*sigma^2)*t)/(sigma*sqrt(t))
}Delta
Call : \({e^{-qt} * N(d_1)}\)
Put : \({e^{-qt} * (N(d_1) - 1 )}\)
getDelta <- function(spot,k,t,rd,rf,sigma,putCallFlag){
d1 <- getD1(spot,k,rd,rf,sigma,t)
if(putCallFlag == "put"){
delta <- -exp(-rf*t)*pnorm(-d1)
}else if(putCallFlag == "call"){
delta <- exp(-rf*t)*pnorm(d1)
}
delta
}Gamma
- Call/Put : \({\frac{e^{-qt} } {(S_0 * \sigma * \sqrt{t})} * N(d_1) }\)
getGamma <- function(spot,k,t,rd,rf,sigma){
d1 <- getD1(spot,k,rd,rf,sigma,t)
exp(-rf*t)*dnorm(d1)/(spot*sigma*sqrt(t))
}\({d_2}\)
\({d_2}\) = \({d_1 - \sigma * \sqrt{t}}\)
Theta \(\theta\)
Call : \(- \frac{S_0 * \sigma * N(d_1)} { 2 * \sqrt{t}} - r * K * {e^{-r * t} * N(d_2)}\)
Put : \(- \frac{S_0 * \sigma * N(d_1)} { 2 * \sqrt{t}} + r * K * {e^{-r * t} * N(-d_2)}\)
getTheta <- function(spot,k,t,rd,rf,sigma,putCallFlag){
d1 <- getD1(spot,k,rd,rf,sigma,t)
d2 <- d1 - sigma*sqrt(t)
# Check why function not multiple by rd, compare with formula?
if(putCallFlag == "put"){
kTerm <- k*exp(-rd*t)*pnorm(-d2)
}else if(putCallFlag == "call"){
kTerm <- -k*exp(-rd*t)*pnorm(d2)
}
-(spot*sigma*dnorm(d1)/(2*sqrt(t))) + kTerm
}計算歷史資料報酬率 \({log(S_1 / S_0) = log(S_1) - log(S_0)}\)
- marketData 歷史資料
- valueDate 評價日
- nDays 每隔幾個資料日期取值計算報酬
- nHistroicalData 距離評價日幾個資料日的資料
- nResampled 依照資料列報酬 隨機抽樣取得報酬
getLogReturn <- function(marketData, valueDate, nDays, nHistroicalData, nResampled) {
# 取得評價日 位置
valeDateLocation <- which(as.Date(marketData$Date) == as.Date(valueDate))
if(length(valeDateLocation)==0){
cat("錯誤：評價日即期匯率不在資料中")
return (0)
}
#取得 評價日價格
spot <- marketData[valeDateLocation,]$Rate
#取得評價日之前的資料列
marketDataByValueDate <- subset(marketData, Date < valueDate)
if(nHistroicalData > length(marketDataByValueDate$Rate)){
cat("錯誤：資料數量不足")
return (0)
}
#取得評價日，指定歷史資料日，每隔幾個資料天 的資料位置
targetData <- seq(valeDateLocation,valeDateLocation + nHistroicalData,nDays)
#取得資料列
spotData <- marketData$Rate[targetData]
#calculate log return, log (s1/s0) = log(s1) - log(s0)
if(!is.null(nResampled)){
if(nResampled < length(spotData)){
cat("警告 : bootstrap 次數 < 選取資料長度")
}
returnRates <- sample(diff(log(spotData)),nResampled,replace = TRUE)
}else{
returnRates <- diff(log(spotData))
}
return (returnRates)
}計算條件風險值，找不到計算公式說明
calculatCVaR <- function(lossVec,extremeValue,confidenceInterval){
nData <- length(lossVec)
alpha <- confidenceInterval/100
fAlpha <- sum(lossVec<=extremeValue)/nData
lambda <- (fAlpha - alpha)/(1-alpha)
CVaRPlus <- ifelse(lambda<1,mean(lossVec[lossVec>extremeValue]),0)
lambda*extremeValue + (1-lambda)*CVaRPlus
}使用RND套件，compute.implied.volatility 計算隱含波動度
if(!require("RND")){
install.packages("RND")
require("RND")
}
setwd("/Users/hadoop/Dropbox/R261進階班/RHandOut/VaRCode/")
#計算信賴度的限值
source("interpolateValue.R")計算VaR
getVarAndEmpiricalDist <- function(marketDate,valueDate,nDays,nHistroicalData,
confidenceInterval,nResampled = NULL,
optionPrice = NULL,putCallFlag = NULL, strike = NULL,
rd = NULL,rf = NULL,tau = NULL, addTheta = FALSE,
dayCounts = 365){
#require("RQuantLib")
require("RND")
source("interpolateValue.R")
# 依照指定評價日與資料列設定，計算log報酬率
returnRates <- getLogReturn (marketData , valueDate , nDays, nHistroicalData, nResampled )
#取得評價日 價格
valueDate
spot <- subset(marketData, as.Date(Date) == as.Date(valueDate))$Rate
spot
str(marketData$Date)
#將log報酬遞減排序
sortR <- sort(returnRates,decreasing = TRUE)
#
nR <- length(returnRates)
# empirical cdf
getPercantage <- function(r,sortPrice){
sum(sortPrice <= r)/nR
}
#依照log 報酬率 計算所有可能的隔日價格
nextSpot <- exp(sortR)*spot
spotloss <- spot - nextSpot
sCdf <- mapply(getPercantage,r = spotloss,MoreArgs = list(sortPrice = spotloss))
extremeValue <- interpolateValue(spotloss,sCdf,confidenceInterval)
#繪製 spot loss 的直方圖
hist(x = spotloss, main = "AUD/USD Spot",
xlab="Spot Loss",ylab = "Frequency", cex.lab= 1.5)
par(new = TRUE)
plot(x = spotloss, y = sCdf, xaxt='n', yaxt='n',type = "S",
col = "blue",xlab = "",ylab="",main="",lwd = 2)
axis(4, col = "blue", col.axis = "blue", lwd = 2)
text(x = spotloss[round(0.02*nR)],y = 0.9*max(sCdf),
labels = "Loss CDF", col = "blue", cex = 1.2)
abline(v = extremeValue ,col = "red", lwd = 2)
##
spotGainLoss <- ifelse(extremeValue>=0,"損失","獲利")
spotCVaR <- calculatCVaR(lossVec = spotloss,
extremeValue = extremeValue,
confidenceInterval = confidenceInterval)
spotCVaRGainLoss <- ifelse(spotCVaR>=0,"損失","獲利")
cat("\n 信賴區間為 ",confidenceInterval,"% 之即期匯率風險值為：",spotGainLoss," ",
abs(extremeValue),sep = "")
cat("\n 即期匯率條件風險值為： ",spotCVaRGainLoss," ",abs(spotCVaR),"\n",sep = "")
# ?p?G???????v?????A?p???ӳ??쭷?I??
if(!is.null(putCallFlag)){
# 計算 Sigma
#sigma <- EuropeanOptionImpliedVolatility(type=putCallFlag,
# value=optionPrice,
# underlying=spot,
# strike=strike,
# dividendYield=rf,
# riskFreeRate=rd,
# maturity=tau,
# volatility=0.1)[1]
sigma <- compute.implied.volatility(r=rd,
te = tau,
s0 = spot,
k = strike,
y=rf,
call.price = optionPrice,
lower = 10^-10,
upper = 0.999)
# 計算 delta
delta <- getDelta(spot = spot,
k = strike,
t = tau,
rd = rd,
rf = rf,
sigma = sigma,
putCallFlag = putCallFlag)
# 計算gamma
gamma <- getGamma(spot = spot,
k = strike,
rd = rd,
rf = rf,
sigma = sigma,
t = tau)
# 計算 theta
if(addTheta){
theta <- getTheta(spot = spot,
k = strike,
t = tau,
rd = rd,
rf = rf,
sigma = sigma,
putCallFlag = putCallFlag)
}else{
theta <- 0
}
# ?p???j?ѷl?q?? empirical distribution
OptionLoss <- -(delta*(nextSpot - spot) + 0.5*gamma*(nextSpot - spot) ^2
+ theta*nDays/dayCounts)
sortV <- sort(OptionLoss)
vCdf <- mapply(getPercantage,r = sortV,MoreArgs = list(sortPrice = sortV))
# ?????v???I??ø??
#windows()
#quartz()
hist(x = sortV, main = paste("Plain Vanilla",putCallFlag),
xlab="Option Loss",ylab = "Frequency", cex.lab= 1.5)
par(new = TRUE)
plot(x = sortV, y = vCdf, xaxt='n', yaxt='n',type = "S",
col = "blue",xlab = "",ylab="",main="",lwd = 2)
axis(4, col = "blue", col.axis = "blue", lwd = 2)
text(x = sortV[round(0.99*nR)],y = 0.9*max(vCdf),
labels = "Loss CDF", col = "blue", cex = 1.2)
vExtremeValue <- interpolateValue(sortV,vCdf,confidenceInterval)
abline(v = vExtremeValue ,col = "red", lwd = 2)
optionGainLoss <- ifelse(vExtremeValue>=0,"損失","獲利")
optionCVaR <- calculatCVaR(lossVec = sortV,
extremeValue = vExtremeValue,
confidenceInterval = confidenceInterval)
#
optionCVaRGainLoss <- ifelse(optionCVaR>=0,"損失","獲利")
cat("\n 信賴區間為 ",confidenceInterval,"% 之選擇權風險值為：",optionGainLoss," ",
abs(vExtremeValue),sep = "")
cat("\n 條件期望風險值為：",optionCVaRGainLoss," ",abs(optionCVaR),sep = "")
cat("\n 選擇權隱含波動度：",sigma,sep = "")
}
}setwd("/Users/hadoop/Dropbox/R261進階班/RHandOut/VaRCode/")
marketData <- read.csv("AUDUSD.csv",header = TRUE, stringsAsFactors = FALSE)
head(marketData)## Date Rate
## 1 2015/8/14 0.734521
## 2 2015/8/13 0.732471
## 3 2015/8/12 0.735541
## 4 2015/8/11 0.740225
## 5 2015/8/10 0.738878
## 6 2015/8/7 0.734817getVarAndEmpiricalDist(marketDate,
valueDate = '2015-08-14',
nDays = 1,
nHistroicalData = 60,
confidenceInterval = 98,
nResampled = 100,
optionPrice = 0.1,
putCallFlag = "call",
strike = 0.7,
rd = 0.004,
rf = 0.0015,
tau = 1.2,
addTheta = TRUE)## chr [1:4075] "2015/8/14" "2015/8/13" "2015/8/12" "2015/8/11" ...
##
## 信賴區間為 98% 之即期匯率風險值為：損失 0.007454164
## 即期匯率條件風險值為： 損失 0.007454164
##
## 信賴區間為 98% 之選擇權風險值為：損失 0.005684087
## 條件期望風險值為：損失 0.005684087
## 選擇權隱含波動度：0.2582577=====================================
f(y+x) = f(y) + f’(y) x + \({\frac{1}{2}}\) f’’(y) \({x^2}\) + \(\epsilon(2)\)
f’(y)：Delta 一階微分
f’’(y)：Gamma 二階微分
\(\epsilon(2)\): 二階近似之誤差