Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- library(reshape2)
- library(ggplot2)
- #Generowanie trajektorii i odczytanie sigmy z danych historycznych
- GBM.simulation<-function(Path, trajectories_number, steps){
- Wig_data<-read.csv(Path)
- Historical_data<- Wig_data[,5]
- Returns <-log(Historical_data[2:length(Historical_data)])-log(Historical_data[1:(length(Historical_data)-1)])
- spot <- Historical_data[length(Historical_data)]
- drift <- mean(Returns)
- sigma <- sd(Returns)
- dt <- 1
- Ncols <- steps
- Nrows <- trajectories_number
- simu.paths <- matrix((drift-sigma^2/2)*dt+sigma*sqrt(dt)*rnorm(Ncols*Nrows), ncol=Ncols, nrow=Nrows)
- simu.paths <- spot*exp( t(apply(simu.paths, 1, cumsum)) )
- x <- matrix(spot,Nrows,1)
- simu.paths<- cbind(x,simu.paths)
- sigma<-sigma*sqrt(252)
- return (list(simu.paths, sigma))
- }
- #Wycena opcji zgodnie z wzorem Blacka-Scholesa
- option.price<-function(spot,strike,sigma,riskfree,tau,type){
- d1<-(log(spot/strike)+(riskfree+sigma^2/2)*tau)/(sigma*sqrt(tau))
- d2<-d1-sigma*sqrt(tau)
- Nd1<-pnorm(d1)
- Nd2<-pnorm(d2)
- discount<-exp(-riskfree*tau)
- price<-spot*Nd1-strike*Nd2*discount
- if(type=="c")
- return(price)
- else if(type=="p")
- return(price-spot+strike*discount)
- }
- #Delta na podstawie wzoru Blacka Scholesa
- delta<-function(spot,strike,sigma,riskfree,tau,type){
- d1<-(log(spot/strike)+(riskfree+sigma^2/2)*tau)/(sigma*sqrt(tau))
- Nd1<-pnorm(d1)
- greek<-Nd1
- if(type=="c")
- return(greek)
- else if(type=="p")
- return(greek-1)
- }
- delta_hedging<-function(simu.paths, trajectories_number, steps, strike, sigma, riskfree, Tau, type, rehedging_freq){
- Nrows<-trajectories_number
- Ncols<-floor(steps/rehedging_freq)
- #Wyznaczenie delty dla kazdego dnia, w ktorym rehedgujemy
- simu.delta<-matrix(0,Nrows,Ncols)
- for(i in 1:Nrows)
- for(k in 1:Ncols)
- simu.delta[i,k]<-delta(spot = simu.paths[i,(1+k*rehedging_freq)], strike, sigma, riskfree, tau = Tau*(1-(k*rehedging_freq)/179), type)
- delta_0<- matrix(delta(spot = simu.paths[1,1], strike, sigma, riskfree, tau = Tau, type),Nrows,1)
- simu.delta<- cbind(delta_0,simu.delta)
- dif.hedge<-apply(simu.delta,1,diff)
- dif.hedge<-t(dif.hedge)
- #Macierz ze stanem gotowki dla kazdego dnia w kazdej z trajektorii
- P.L<-matrix(0,Nrows,(Ncols+1))
- for(i in 1:Nrows)
- 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]
- for(j in 1:Nrows)
- for(k in 2:(Ncols+1))
- 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)]
- #Gotowka w dniu zapadniecia opcji
- final.position<-P.L[,(Ncols+1)]*exp(((steps-1)%%rehedging_freq)*riskfree/252)
- #Wartosc WIG20 w dniu zapadniecia opcji
- final.price<-simu.paths[,steps]
- #Koncowy zysk/strata (mamy gotowke, zamykamy krotka/dluga pozycje na indeksie, opcja wykonuje sie lub nie)
- final.pl <- rep(0,Nrows)
- for(i in 1:Nrows){
- if(type == "p" & final.price[i] < strike)
- final.pl[i] <- final.position[i] + (simu.delta[i,1] + sum(dif.hedge[i,]))*simu.paths[i,steps] - (strike - simu.paths[i,steps])
- else if(type =="c" & final.price[i] > strike)
- final.pl[i] <- final.position[i] + (simu.delta[i,1] + sum(dif.hedge[i,]))*simu.paths[i,steps] - (simu.paths[i,steps] - strike)
- else {
- final.pl[i] <- final.position[i] + (simu.delta[i,1] + sum(dif.hedge[i,]))*simu.paths[i,steps]
- }
- }
- #zakladam rehedging_freq rozne od 1
- cash<-matrix(0,Nrows,(steps+1))
- WIG20<-matrix(0,Nrows,(steps+1))
- if(rehedging_freq!=1){
- for(j in 1:(Ncols+1))
- cash[,(1+(j-1)*rehedging_freq)]<-P.L[,j]
- for(k in 1:Ncols)
- for(l in 1:(rehedging_freq-1))
- cash[,(1+(k-1)*rehedging_freq+l)]<-cash[,(1+(k-1)*rehedging_freq)]*exp(l*riskfree/252)
- m<-(steps-1)%%rehedging_freq #moze byc 0
- for(n in 1:m)
- cash[,(1+ Ncols*rehedging_freq+n)]<-cash[,(1+Ncols*rehedging_freq)]*exp(n*riskfree/252)
- cash[,(steps+1)]<-final.pl
- #indeks WIG20
- for(o in 1:(Ncols+1))
- WIG20[,(1+(o-1)*rehedging_freq)]<-simu.delta[,o]*simu.paths[,(1+(o-1)*rehedging_freq)]
- for(p in 1:Ncols)
- for(q in 1:(rehedging_freq-1))
- WIG20[,(1+(p-1)*rehedging_freq+q)]<-simu.delta[,p]*simu.paths[,(1+(p-1)*rehedging_freq+q)]
- for(r in 1:m)
- WIG20[,(1+Ncols*rehedging_freq+r)]<-simu.delta[,(Ncols+1)]*simu.paths[,(1+ Ncols*rehedging_freq+r)]
- WIG20[,(steps+1)]<-0
- }
- else{
- cash<-P.L
- cash[,(steps+1)]<-final.pl
- for(s in 1:Nrows)
- WIG20[s,]<-simu.delta[s,]*simu.paths[s,]
- WIG20[,(steps+1)]<-0
- }
- return (list(final.pl, cash, WIG20))
- }
- #Parametr 1 - liczba symulowanych dni
- steps <- 179
- #Parametr 2 - liczba generowanych trajektorii
- trajectories_number <- 1000
- x<-GBM.simulation("wig20.csv", 10 * trajectories_number, steps)
- #Parametr 3 - stopa wolna od ryzyka
- riskfree <- 0.0515
- #Parametr 4 - strike opcji
- strike <- c(1900, 2600, 3200)
- #Parametr 5 - czas do wygasniecia (dla opcji zapadajacych 16 wrzesnia to bedzie 179 dni)
- tau <- 179/252
- #Parametr 6 - rodzaj opcji "p" - put "c" - call
- type <- c("p", "c")
- #Parametr 7 - częstotliwość aktualizowania portfela (1 - codziennie, 5 - co tydzien, 21 - co miesiac)
- rehedging_freq <- c(1, 5, 10, 21)
- #Koncowa wartosc portfela dla kazdej wygenerowanej trajektorii
- #profit_loss<- delta_hedging(simu.paths = x[[1]],trajectories_number, steps, strike[3], sigma = x[[2]], riskfree, tau, type[2], rehedging_freq[2])
- #p_l<-profit_loss[[1]]
- #cash<-profit_loss[[2]]
- #index<-profit_loss[[3]]
- profit_loss_matrix <- matrix(nrow = 1000, ncol = 4)
- profit_loss_source <- data.frame()
- chart <- list(ggplot())
- for(k in 1:6) {
- for(i in 1:4) {
- profit_loss <- delta_hedging(simu.paths = x[[1]], trajectories_number,
- steps, strike[((k - 1) %% 3) + 1], sigma = x[[2]], riskfree,
- tau, type[floor(1 + (k - 1) / 3)], rehedging_freq[i])[[1]]
- profit_loss_matrix[ , i] <- t(profit_loss)
- }
- profit_loss_source <- ggplot.source(profit_loss_matrix, y_values_name = 'profit',
- types = c('codziennie', 'co tydzien', 'co dwa tygodnie', 'co miesiac'),
- types_name = 'rehedging_frequency')
- if(type[floor(1 + (k - 1) / 3)] == 'p')
- title <- 'put'
- else if(type[floor(1 + (k - 1) / 3)] == 'c')
- title <- 'call'
- binnr <- 1 + 3.3 * log(trajectories_number)
- binwh <- (max(profit_loss_matrix) - min(profit_loss_matrix)) / binnr
- chart[[k]] <- ggplot(profit_loss_source, aes(x = profit, fill = rehedging_frequency)) +
- geom_histogram(aes(colour = rehedging_frequency), position = 'identity', bins = binnr, binwidth = binwh) +
- facet_wrap(~rehedging_frequency, nrow = 2) + theme_light() +
- ggtitle(paste(title, '@', strike[((k - 1) %% 3) + 1], sep = ""))
- ggsave(paste('D:/Studia/Inżynieria finansowa 1/Projekt/Wykresy/rehedging_',
- title,
- strike[((k - 1) %% 3) + 1],
- '.png',
- sep = ''),
- dpi = 1800)
- }
- profit_q <- matrix(nrow = 10001, ncol = 4)
- for(k in 1:6) {
- for(i in 1:4) {
- profit_loss <- delta_hedging(simu.paths = x[[1]], 10 * trajectories_number,
- steps, strike[((k - 1) %% 3) + 1], sigma = x[[2]], riskfree,
- tau, type[floor(1 + (k - 1) / 3)], rehedging_freq[i])[[1]]
- profit_q[, i] <- quantile(profit_loss, probs = seq(0, 1, by = 0.0001))
- }
- profit_q_source <- ggplot.source(profit_q, x_values = seq(0, 1, by = 0.0001),
- x_values_name = 'quantile', y_values_name = 'profit',
- types = c('codziennie', 'co tydzien', 'co dwa tygodnie', 'co miesiac'),
- types_name = 'rehedging_frequency')
- if(type[floor(1 + (k - 1) / 3)] == 'p')
- title <- 'put'
- else if(type[floor(1 + (k - 1) / 3)] == 'c')
- title <- 'call'
- chart[[k]] <- ggplot(profit_q_source, aes(x = quantile, y = profit, fill = rehedging_frequency)) +
- geom_line(aes(colour = rehedging_frequency)) +
- theme_light() +
- ggtitle(paste(title, '@', strike[((k - 1) %% 3) + 1], sep = ""))
- ggsave(paste('D:/Studia/Inżynieria finansowa 1/Projekt/Wykresy/rehedging_quantiles',
- title,
- strike[((k - 1) %% 3) + 1],
- '.png',
- sep = ''),
- dpi = 1800)
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement