API tools faq
paste
Guest User

Untitled

a guest
Nov 22nd, 2017
81
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 8.67 KB | None | 0 0
raw download clone embed print report
  1. # Load all functions related to the experiement
  2. source("functions_simulation_results.R")
  3.  
  4. #load simulation parameters
  5. source("simulation_parameters.R")
  6.  
  7. # Generated 500 1-day paths
  8. source("model_parameters.R")
  9. source("observation_simulation_euler_P.R")
  10.  
  11. # SELECT WITH AND WITHOUT JUMPS + SAMPLE SIZE
  12. S$jump.v <- FALSE
  13. S$jump.y <- FALSE
  14. l.sample <- 500
  15.  
  16. # get indices of paths with either a jump in V or Y
  17. ijyv <- unique(which(Z.y != 0 | Z.v != 0, arr.ind = TRUE)[,1])
  18. # get indices of paths without any jumps
  19. inoj <- as.vector(1:S$n.traj.obs)[-ijyv]
  20.  
  21. # sample 500 paths with/without jumps
  22. if(S$jump.v == FALSE | S$jump.y == FALSE){
  23. i <- sample(x = inoj,size = l.sample,replace = FALSE)
  24. V <- V[i,]
  25. Y <- Y[i,]
  26. Z.y <- Z.y[i,]
  27. Z.v <- Z.v[i,]
  28. dI <- dI[i,]
  29. }else{
  30. i <- sample(x = ijyv,size = l.sample,replace = FALSE)
  31. V <- V[i,]
  32. Y <- Y[i,]
  33. Z.y <- Z.y[i,]
  34. Z.v <- Z.v[i,]
  35. dI <- dI[i,]
  36. }
  37.  
  38.  
  39. # option related parameters
  40. BSdelta <- S$delta
  41. is.call <- 1
  42. tau <- 30/252
  43. M <- c(1,2,3,5,10,1560)
  44.  
  45. lm <- length(M)
  46. rmse.orv <- matrix(ncol = S$ld, nrow = lm)
  47. rrmse.orv <- matrix(ncol = S$ld, nrow = lm)
  48. r2.orv <- matrix(ncol = S$ld, nrow = lm)
  49. rmse.ivrv <- matrix(ncol = S$ld, nrow = lm)
  50. rrmse.ivrv <- matrix(ncol = S$ld, nrow = lm)
  51. r2.ivrv <- matrix(ncol = S$ld, nrow = lm)
  52.  
  53. for(o in 1:S$ld){
  54.  
  55. tito <- as.character(o)
  56.  
  57. # Get end of day strikes for each path
  58. GetMoneyness <- function(V, t, r, q, is.call, BSdelta){
  59. fun <- function(x) GetDelta(x, V=V, t=t,
  60. r=r, q=q, is.call = is.call) - BSdelta
  61. moneyness <- uniroot(f = fun, interval = c(0,2))$root
  62. return(moneyness)
  63. }
  64.  
  65. moneyness <- mapply(FUN=GetMoneyness, V=V[,ncol(V)],
  66. MoreArgs = list(t=tau,
  67. r=P$r,
  68. q=P$q,
  69. is.call=is.call,
  70. BSdelta=BSdelta[o]))
  71.  
  72. K <- exp(Y[,ncol(Y)])/moneyness
  73. K.mat <- matrix(rep(K,ncol(Y)),nrow = nrow(Y))
  74.  
  75. # info for grid interpolation
  76. S0 <- exp(Y)
  77. moneyness.intraday <- S0/K
  78. V0 <- V
  79. logV <- log(V0)
  80. loc <- cbind( as.vector(logV), as.vector(moneyness.intraday))
  81.  
  82. # Get IV and price through grid interpolation
  83. IV <- matrix(interp.surface( S$obj.30.iv, loc), ncol=S$n.step.obs+1)
  84. prices <- GetBlackScholesCall(S=S0,
  85. K=K.mat ,
  86. t=tau,
  87. r=P$r,
  88. q=P$q,
  89. sigma=IV,
  90. is.call=is.call)
  91.  
  92. # Get Greeks through grid interpolation
  93. delta <- matrix(interp.surface( S$obj.30.d, loc), ncol=S$n.step.obs+1)*S0
  94. vega <- matrix(interp.surface( S$obj.30.v, loc), ncol=S$n.step.obs+1)*S0
  95.  
  96. # compute IVQV derivative
  97. d <- (log(S0/K) + (P$r - P$q + 0.5*IV^2) * tau)/
  98. (IV * sqrt(tau))
  99. phi.d <- dnorm(d, 0, 1)
  100. inverse.BS.vega <- 1/(S0 * phi.d * sqrt(tau))
  101. dIV.y <- delta * inverse.BS.vega
  102. dIV.v <- vega * inverse.BS.vega
  103.  
  104. # Get cumulative sum values
  105. Z.y.cs <- t(apply(Z.y, 1, cumsum))
  106. Z.v.cs <- t(apply(Z.v, 1, cumsum))
  107. dI.cs <- t(apply(dI, 1, cumsum))
  108.  
  109. #For each M and each option delta, compute OQV/IVQV
  110. for(m in 1:lm){
  111.  
  112. titm <- as.character(M[m])
  113. title <- paste0("M = ", titm)
  114.  
  115. # depending on M used, get beginning of period indices and end of period indices
  116. if(m == 1){
  117. ide <- 1560
  118. idb <- 1
  119. }else{
  120. ide <- seq(from=S$n.step.obs/M[m], to = S$n.step.obs, by = S$n.step.obs/M[m])
  121. idb <- seq(1, to = S$n.step.obs, by = S$n.step.obs/M[m] )
  122. }
  123.  
  124. Z.y.mat <- cbind(rep(0, l.sample), Z.y.cs[,ide])
  125. Z.v.mat <- cbind(rep(0, l.sample), Z.v.cs[,ide])
  126. dI.mat <- cbind(rep(0, l.sample), dI.cs[,ide])
  127. dZ.y.mat <- matrix(apply(Z.y.mat, 1, diff),nrow=l.sample, byrow=TRUE)
  128. dZ.v.mat <- matrix(apply(Z.v.mat, 1, diff),nrow=l.sample, byrow=TRUE)
  129. I.mat <- matrix(apply(dI.mat, 1, diff),nrow=l.sample, byrow=TRUE)
  130.  
  131. # Get before jumps values for interpolation
  132. S0_ <- exp(Y[,ide+1]-dZ.y.mat)
  133. moneyness.intraday_ <- S0_/K.mat[,ide+1]
  134. V0_ <- V[,ide+1]-dZ.v.mat
  135. V0_[which(V0_ < P$vmin)] <- P$vmin
  136. logV_ <- log(V0_)
  137. loc_ <- cbind( as.vector(logV_),
  138. as.vector(moneyness.intraday_))
  139.  
  140. # Get IV and prices before jumps
  141. IV_ <- matrix(interp.surface( S$obj.30.iv, loc_), ncol=ncol(S0_))
  142. prices_ <- GetBlackScholesCall(S=S0_,
  143. K=K.mat[,ide+1] ,
  144. t=tau,
  145. r=P$r,
  146. q=P$q,
  147. sigma=IV_,
  148. is.call=is.call)
  149.  
  150. #OQV
  151. mat1 <- (delta[,idb]^2 + (P$sigma * vega[,idb])^2 +
  152. 2*P$sigma*P$rho*delta[,idb]*vega[,idb])* I.mat
  153. term1 <- apply(mat1 ,1,sum)
  154. mat2 <- matrix((prices[,ide+1] - prices_)^2, nrow=l.sample)
  155. term2 <- apply(mat2, 1, sum)
  156. OQV <- (term1 + term2)
  157. ORV <- apply(prices, 1, GetRVSimulation)
  158.  
  159. # Scatterplot of ORV/OQV
  160. max <- max(OQV)
  161. ran <- c(0,max)
  162. plot(x = ORV, y = OQV,xlab = "ORV", ylab = "OQV", main= title,
  163. xlim=ran, ylim=ran,pch=19)
  164. lines(x=ran,y=ran)
  165. imgname <- paste0("scatterplot_ORV_",tito,"_", titm, ".png")
  166. dev.copy(png, imgname)
  167. dev.off()
  168.  
  169. # Regression ORV/OQV
  170. rmse.orv[m,o] <- sqrt (mean((OQV-ORV)^2))
  171. rrmse.orv[m,o] <- (sqrt(sum(OQV-ORV)^2))/
  172. (length(ORV) * mean(ORV))
  173. oqv.fit <- lm(OQV ~ ORV)
  174. r2.orv[m,o] <- summary(oqv.fit)$r.squared
  175.  
  176. # IVQV
  177. mat1 <- (dIV.y[,idb]^2 + (P$sigma * dIV.v[,idb])^2 +
  178. 2*P$sigma*P$rho*dIV.y[,idb]*dIV.v[,idb])* I.mat
  179. term1 <- apply(mat1 ,1,sum)
  180. mat2 <- matrix((IV[,ide+1] - IV_)^2, nrow=l.sample)
  181. term2 <- apply(mat2, 1, sum)
  182.  
  183. IVQV <- (term1 + term2 )
  184. idx.ignore <- which(is.nan(IVQV) | IVQV > 0.1) # when the option is too deep out of the money, inverse BS vega = Inf, and IVQV is NaN
  185.  
  186. if(length(idx.ignore) > 0){
  187. IVQV <- IVQV[-idx.ignore]
  188. IVRV <- apply(IV, 1, GetRVSimulation)[-idx.ignore]
  189. }else{
  190. IVQV <- IVQV
  191. IVRV <- apply(IV, 1, GetRVSimulation)
  192. }
  193.  
  194. # Scatterplot of IVRV/IVQV
  195. max <- max(IVQV)
  196. ran <- c(0,max)
  197. plot(x = IVRV, y = IVQV,xlab = "IVRV", ylab = "IVQV", main= title,
  198. xlim=ran, ylim=ran,pch=19)
  199. lines(x=ran,y=ran)
  200. titm <- as.character(M[m])
  201. imgname <- paste0("scatterplot_IVRV_",tito,"_", titm, ".png")
  202. dev.copy(png, imgname)
  203. dev.off()
  204.  
  205. # Regression IVRV/VQV
  206. rmse.ivrv[m,o] <- sqrt (mean((IVQV-IVRV)^2))
  207. rrmse.ivrv[m,o] <- (sqrt(sum(IVQV-IVRV)^2))/
  208. (length(IVRV) * mean(IVRV))
  209. ivqv.fit <- lm(IVQV ~ IVRV)
  210. r2.ivrv[m,o] <- summary(ivqv.fit)$r.squared
  211.  
  212. }
  213.  
  214.  
  215.  
  216. }
  217.  
  218. colname <- as.character(S$delta)
  219. rowname <- as.character(M)
  220.  
  221. colnames(rmse.ivrv) <- colname
  222. colnames(rmse.orv) <- colname
  223. colnames(rrmse.ivrv)<- colname
  224. colnames(rrmse.orv) <- colname
  225. colnames(r2.ivrv) <- colname
  226. colnames(r2.orv) <- colname
  227.  
  228. row.names(rmse.ivrv) <- rowname
  229. row.names(rmse.orv) <- rowname
  230. row.names(rrmse.ivrv)<- rowname
  231. row.names(rrmse.orv) <- rowname
  232. row.names(r2.ivrv) <- rowname
  233. row.names(r2.orv) <- rowname
  234.  
  235. results.ivrv <- rbind(rmse.ivrv,rrmse.ivrv,r2.ivrv)
  236. results.orv <- rbind(rmse.orv,rrmse.orv,r2.orv)
  237.  
  238. write.table(x= results.ivrv, file = "Results_IVRV_withoutjumps.csv", sep = ",", col.names = colname)
  239. write.table(x= results.orv, file = "Results_ORV_withoutjumps.csv", sep = ",", col.names = colname)
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!