API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Dec 3rd, 2016
78
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
R 8.70 KB | None | 0 0
raw download clone embed print report
  1. library(reshape2)
  2. library(ggplot2)
  3.  
  4. #Generowanie trajektorii i odczytanie sigmy z danych historycznych
  5.  
  6. GBM.simulation<-function(Path, trajectories_number, steps){
  7.  
  8.   Wig_data<-read.csv(Path)
  9.   Historical_data<- Wig_data[,5]
  10.  
  11.   Returns <-log(Historical_data[2:length(Historical_data)])-log(Historical_data[1:(length(Historical_data)-1)])
  12.   spot <- Historical_data[length(Historical_data)]
  13.   drift <- mean(Returns)
  14.   sigma <- sd(Returns)
  15.   dt <- 1
  16.  
  17.   Ncols <- steps
  18.   Nrows <- trajectories_number
  19.  
  20.   simu.paths <- matrix((drift-sigma^2/2)*dt+sigma*sqrt(dt)*rnorm(Ncols*Nrows), ncol=Ncols, nrow=Nrows)
  21.   simu.paths <- spot*exp( t(apply(simu.paths, 1, cumsum)) )
  22.   x <- matrix(spot,Nrows,1)
  23.   simu.paths<- cbind(x,simu.paths)
  24.  
  25.   sigma<-sigma*sqrt(252)
  26.  
  27.   return (list(simu.paths, sigma))
  28. }
  29.  
  30. #Wycena opcji zgodnie z wzorem Blacka-Scholesa
  31.  
  32. option.price<-function(spot,strike,sigma,riskfree,tau,type){
  33.  
  34.   d1<-(log(spot/strike)+(riskfree+sigma^2/2)*tau)/(sigma*sqrt(tau))
  35.   d2<-d1-sigma*sqrt(tau)
  36.   Nd1<-pnorm(d1)
  37.   Nd2<-pnorm(d2)
  38.   discount<-exp(-riskfree*tau)
  39.   price<-spot*Nd1-strike*Nd2*discount
  40.  
  41.   if(type=="c")
  42.     return(price)
  43.   else if(type=="p")
  44.     return(price-spot+strike*discount)
  45. }
  46.  
  47. #Delta na podstawie wzoru Blacka Scholesa
  48.  
  49. delta<-function(spot,strike,sigma,riskfree,tau,type){
  50.  
  51.   d1<-(log(spot/strike)+(riskfree+sigma^2/2)*tau)/(sigma*sqrt(tau))
  52.   Nd1<-pnorm(d1)
  53.   greek<-Nd1
  54.  
  55.   if(type=="c")
  56.     return(greek)
  57.   else if(type=="p")
  58.     return(greek-1)
  59. }
  60.  
  61.  
  62. delta_hedging<-function(simu.paths, trajectories_number, steps, strike, sigma, riskfree, Tau, type, rehedging_freq){
  63.  
  64.   Nrows<-trajectories_number
  65.   Ncols<-floor(steps/rehedging_freq)
  66.  
  67.   #Wyznaczenie delty dla kazdego dnia, w ktorym rehedgujemy
  68.  
  69.   simu.delta<-matrix(0,Nrows,Ncols)
  70.   for(i in 1:Nrows)
  71.     for(k in 1:Ncols)
  72.       simu.delta[i,k]<-delta(spot = simu.paths[i,(1+k*rehedging_freq)], strike, sigma, riskfree, tau = Tau*(1-(k*rehedging_freq)/179), type)
  73.  
  74.   delta_0<- matrix(delta(spot = simu.paths[1,1], strike, sigma, riskfree, tau = Tau, type),Nrows,1)
  75.   simu.delta<- cbind(delta_0,simu.delta)
  76.  
  77.   dif.hedge<-apply(simu.delta,1,diff)
  78.   dif.hedge<-t(dif.hedge)
  79.  
  80.  
  81.   #Macierz ze stanem gotowki dla kazdego dnia w kazdej z trajektorii
  82.  
  83.   P.L<-matrix(0,Nrows,(Ncols+1))
  84.  
  85.   for(i in 1:Nrows)
  86.     P.L[i,1]<- option.price(spot = simu.paths[i,1], strike, sigma, riskfree, tau = Tau, type) - simu.delta[i,1]*simu.paths[i,1]
  87.  
  88.   for(j in 1:Nrows)
  89.     for(k in 2:(Ncols+1))
  90.       P.L[j,k]<-P.L[j,k-1]*exp(riskfree*rehedging_freq/252) - dif.hedge[j,k-1]*simu.paths[j,(1+(k-1)*rehedging_freq)]
  91.  
  92.   #Gotowka w dniu zapadniecia opcji
  93.   final.position<-P.L[,(Ncols+1)]*exp(((steps-1)%%rehedging_freq)*riskfree/252)
  94.  
  95.   #Wartosc WIG20 w dniu zapadniecia opcji
  96.   final.price<-simu.paths[,steps]
  97.  
  98.   #Koncowy zysk/strata (mamy gotowke, zamykamy krotka/dluga pozycje na indeksie, opcja wykonuje sie lub nie)
  99.   final.pl <- rep(0,Nrows)
  100.  
  101.   for(i in 1:Nrows){
  102.     if(type == "p" & final.price[i] < strike)
  103.       final.pl[i] <- final.position[i] + (simu.delta[i,1] + sum(dif.hedge[i,]))*simu.paths[i,steps] - (strike - simu.paths[i,steps])
  104.     else if(type =="c" & final.price[i] > strike)
  105.       final.pl[i] <- final.position[i] + (simu.delta[i,1] + sum(dif.hedge[i,]))*simu.paths[i,steps] - (simu.paths[i,steps] - strike)
  106.     else {
  107.       final.pl[i] <- final.position[i] + (simu.delta[i,1] + sum(dif.hedge[i,]))*simu.paths[i,steps]
  108.     }
  109.   }
  110.  
  111.   #zakladam rehedging_freq rozne od 1
  112.  
  113.   cash<-matrix(0,Nrows,(steps+1))
  114.   WIG20<-matrix(0,Nrows,(steps+1))
  115.  
  116.   if(rehedging_freq!=1){
  117.    
  118.     for(j in 1:(Ncols+1))
  119.       cash[,(1+(j-1)*rehedging_freq)]<-P.L[,j]
  120.    
  121.     for(k in 1:Ncols)
  122.       for(l in 1:(rehedging_freq-1))
  123.         cash[,(1+(k-1)*rehedging_freq+l)]<-cash[,(1+(k-1)*rehedging_freq)]*exp(l*riskfree/252)
  124.      
  125.       m<-(steps-1)%%rehedging_freq #moze byc 0
  126.      
  127.       for(n in 1:m)
  128.         cash[,(1+ Ncols*rehedging_freq+n)]<-cash[,(1+Ncols*rehedging_freq)]*exp(n*riskfree/252)
  129.      
  130.       cash[,(steps+1)]<-final.pl
  131.      
  132.       #indeks WIG20
  133.      
  134.      
  135.       for(o in 1:(Ncols+1))
  136.         WIG20[,(1+(o-1)*rehedging_freq)]<-simu.delta[,o]*simu.paths[,(1+(o-1)*rehedging_freq)]
  137.      
  138.       for(p in 1:Ncols)
  139.         for(q in 1:(rehedging_freq-1))
  140.           WIG20[,(1+(p-1)*rehedging_freq+q)]<-simu.delta[,p]*simu.paths[,(1+(p-1)*rehedging_freq+q)]
  141.      
  142.      
  143.       for(r in 1:m)
  144.         WIG20[,(1+Ncols*rehedging_freq+r)]<-simu.delta[,(Ncols+1)]*simu.paths[,(1+ Ncols*rehedging_freq+r)]
  145.      
  146.       WIG20[,(steps+1)]<-0
  147.   }
  148.   else{
  149.    
  150.     cash<-P.L
  151.     cash[,(steps+1)]<-final.pl
  152.    
  153.     for(s in 1:Nrows)
  154.       WIG20[s,]<-simu.delta[s,]*simu.paths[s,]
  155.    
  156.     WIG20[,(steps+1)]<-0
  157.    
  158.   }
  159.   return (list(final.pl, cash, WIG20))
  160. }
  161.  
  162. #Parametr 1 - liczba symulowanych dni
  163. steps <- 179
  164.  
  165. #Parametr 2 - liczba generowanych trajektorii
  166. trajectories_number <- 1000
  167.  
  168. x<-GBM.simulation("wig20.csv", 10 * trajectories_number, steps)
  169.  
  170. #Parametr 3 - stopa wolna od ryzyka
  171. riskfree <- 0.0515
  172.  
  173. #Parametr 4 - strike opcji
  174. strike <- c(1900, 2600, 3200)
  175.  
  176. #Parametr 5 - czas do wygasniecia (dla opcji zapadajacych 16 wrzesnia to bedzie 179 dni)
  177. tau <- 179/252
  178.  
  179. #Parametr 6 - rodzaj opcji "p" - put "c" - call
  180. type <- c("p", "c")
  181.  
  182. #Parametr 7 - częstotliwość aktualizowania portfela (1 - codziennie, 5 - co tydzien, 21 - co miesiac)
  183. rehedging_freq <- c(1, 5, 10, 21)
  184.  
  185. #Koncowa wartosc portfela dla kazdej wygenerowanej trajektorii
  186. #profit_loss<- delta_hedging(simu.paths = x[[1]],trajectories_number, steps, strike[3], sigma = x[[2]], riskfree, tau, type[2], rehedging_freq[2])
  187.  
  188. #p_l<-profit_loss[[1]]
  189. #cash<-profit_loss[[2]]
  190. #index<-profit_loss[[3]]
  191.  
  192. profit_loss_matrix <- matrix(nrow = 1000, ncol = 4)
  193. profit_loss_source <- data.frame()
  194. chart <- list(ggplot())
  195.  
  196. for(k in 1:6) {
  197.  
  198.   for(i in 1:4) {
  199.    
  200.     profit_loss <- delta_hedging(simu.paths = x[[1]], trajectories_number,
  201.                                  steps, strike[((k - 1) %% 3) + 1], sigma = x[[2]], riskfree,
  202.                                  tau, type[floor(1 + (k - 1) / 3)], rehedging_freq[i])[[1]]
  203.    
  204.    
  205.     profit_loss_matrix[ , i] <- t(profit_loss)
  206.    
  207.   }
  208.  
  209.   profit_loss_source <- ggplot.source(profit_loss_matrix, y_values_name = 'profit',
  210.                                       types = c('codziennie', 'co tydzien', 'co dwa tygodnie', 'co miesiac'),
  211.                                       types_name = 'rehedging_frequency')
  212.  
  213.   if(type[floor(1 + (k - 1) / 3)] == 'p')
  214.     title <- 'put'
  215.   else if(type[floor(1 + (k - 1) / 3)] == 'c')
  216.     title <- 'call'
  217.  
  218.   binnr <- 1 + 3.3 * log(trajectories_number)
  219.   binwh <- (max(profit_loss_matrix) - min(profit_loss_matrix)) / binnr
  220.  
  221.   chart[[k]] <- ggplot(profit_loss_source, aes(x = profit, fill = rehedging_frequency)) +
  222.     geom_histogram(aes(colour = rehedging_frequency), position = 'identity', bins = binnr, binwidth = binwh) +
  223.     facet_wrap(~rehedging_frequency, nrow = 2) + theme_light() +
  224.     ggtitle(paste(title, '@', strike[((k - 1) %% 3) + 1], sep = ""))
  225.  
  226.   ggsave(paste('D:/Studia/Inżynieria finansowa 1/Projekt/Wykresy/rehedging_',
  227.                title,
  228.                strike[((k - 1) %% 3) + 1],
  229.                '.png',
  230.                sep = ''),
  231.          dpi = 1800)
  232.  
  233. }
  234.  
  235. profit_q <- matrix(nrow = 10001, ncol = 4)
  236.  
  237. for(k in 1:6) {
  238.  
  239.   for(i in 1:4) {
  240.    
  241.     profit_loss <- delta_hedging(simu.paths = x[[1]], 10 * trajectories_number,
  242.                                  steps, strike[((k - 1) %% 3) + 1], sigma = x[[2]], riskfree,
  243.                                  tau, type[floor(1 + (k - 1) / 3)], rehedging_freq[i])[[1]]
  244.    
  245.    
  246.     profit_q[, i] <- quantile(profit_loss, probs = seq(0, 1, by = 0.0001))
  247.    
  248.   }
  249.  
  250.   profit_q_source <- ggplot.source(profit_q, x_values = seq(0, 1, by = 0.0001),
  251.                                    x_values_name = 'quantile', y_values_name = 'profit',
  252.                                       types = c('codziennie', 'co tydzien', 'co dwa tygodnie', 'co miesiac'),
  253.                                       types_name = 'rehedging_frequency')
  254.  
  255.   if(type[floor(1 + (k - 1) / 3)] == 'p')
  256.     title <- 'put'
  257.   else if(type[floor(1 + (k - 1) / 3)] == 'c')
  258.     title <- 'call'
  259.  
  260.   chart[[k]] <- ggplot(profit_q_source, aes(x = quantile, y = profit, fill = rehedging_frequency)) +
  261.     geom_line(aes(colour = rehedging_frequency)) +
  262.     theme_light() +
  263.     ggtitle(paste(title, '@', strike[((k - 1) %% 3) + 1], sep = ""))
  264.  
  265.   ggsave(paste('D:/Studia/Inżynieria finansowa 1/Projekt/Wykresy/rehedging_quantiles',
  266.                title,
  267.                strike[((k - 1) %% 3) + 1],
  268.                '.png',
  269.                sep = ''),
  270.          dpi = 1800)
  271.  
  272. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
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!