performance_with_gearing_osebx.r

# This source code is provided under the following Creative Commons license: 
#   Attribution-NonCommercial-ShareAlike 4.0 International 
#   https://creativecommons.org/licenses/by-nc-sa/4.0/
# Please contact the author for commercial use: andre@andersen.im

# Fixed parameters

exchange <- 'OSE'
instrument <- 'OSEBX'
p_level = 0.90
allowance <- 20000
sim_count <- 10000
risk_free_return <- 0.0056 # STATSOBL-3Y-EFF Last updated: Thursday 25. August 09:10

# Setting up the different scenarios

destfile <- 'scenarios.csv'
scenarios <- read.csv(file=destfile, strip.white=TRUE)
scenarios$sample_start <- as.Date(scenarios$sample_start)
scenarios$sample_end <- as.Date(scenarios$sample_end)
scenarios$invest_start <- as.Date(scenarios$invest_start) 
scenarios$invest_end <- as.Date(scenarios$invest_end)

# Defining functions for simulation and data processing

build_period <- function(sample_start_date, sample_end_date){
  library(lubridate)
  period <- as.data.frame(seq(sample_start_date, sample_end_date, by="days"))
  names(period) <- 'quote_date'
  period$eom <- period$quote_date == period$quote_date - days(day(period$quote_date)) + 1 + months(1) - 1
  period$som <- period$quote_date == period$quote_date - days(day(period$quote_date)) + 1
  return(period)
}

fetch_dataset <- function(start_date, end_date, exchange, instrument){
  period <- build_period(start_date, end_date)
  url <- paste0('http://www.netfonds.no/quotes/paperhistory.php?paper=',instrument,'.',exchange,'&csv_format=csv')
  destfile <- 'osebx.csv'
  
  if(!file.exists(destfile)){
    download.file(url,destfile)
  }
  
  dataset <- read.csv(file=destfile)
  dataset$quote_date <- as.Date(as.character(dataset$quote_date), "%Y%m%d")
  
  dataset <- dataset[order(dataset$quote_date),]
  
  previous_row_date <- max(dataset$quote_date[dataset$quote_date < start_date])
  
  data_fields <- c('paper','exch','open', 'high', 'low', 'close', 'volume', 'value')
  
  previous_row <- dataset[dataset$quote_date==previous_row_date,data_fields]
  
  dataset <- merge(x=period, y=dataset, by='quote_date', all.x=TRUE)
  
  for(i in 1:nrow(dataset)) {
    row <- dataset[i,]
    if(is.na(row$close) || is.na(row$low)){
      dataset[i,data_fields] <- previous_row
    } else {
      previous_row <- dataset[i, data_fields]
    }
  }
  
  dataset_size <- nrow(dataset)
  ror <- diff(dataset$close)/dataset$close[-dataset_size]
  ror <- append(ror, 0, 0)
  dataset$ror <- ror
  ror_low <- diff(dataset$low)/dataset$low[-dataset_size]
  ror_low <- append(ror_low, 0, 0)
  dataset$ror_low <- ror_low
  dataset <- subset(dataset, select = c('quote_date', 'som', 'eom', 'ror', 'ror_low'))
  return(dataset)
}

sample_df = function(df,n){
  return(df[sample(nrow(df),n),])
}

resample <- function(dataset_invest, dataset_sample){
  sample_df = function(df,n){
    # http://stackoverflow.com/a/8273414/604048
    return(df[sample(nrow(df),n,replace=TRUE),])
  }
  dataset_resampled <- sample_df(dataset_sample, nrow(dataset_invest))
  dataset_resampled$quote_date <- dataset_invest$quote_date
  dataset_resampled$som <- dataset_invest$som
  dataset_resampled$eom <- dataset_invest$eom
  rownames(dataset_resampled) <- 1:nrow(dataset_resampled)
  return(dataset_resampled)
}

simulate <- function(dataset, allowance, gearing_factor, interest_rate){
  loan_flow <- allowance * gearing_factor
  cash_flow <- allowance + loan_flow
  interest_rate_monthly <- ((1+interest_rate) ^ (1/12))-1
  
  dataset$loan_balance <- 0
  dataset$accrued_interest <- 0
  dataset$accrued_allowance <- 0
  dataset$asset_balance <- 0
  dataset$margin_call <- FALSE
  
  loan_balance <- 0
  accrued_interest <- 0
  accrued_allowance <- 0
  asset_balance <- 0
  
  for(i in 1:nrow(dataset)) {
    
    if(dataset[i,'som']){
      asset_balance <- asset_balance + cash_flow
      loan_balance <- loan_balance + loan_flow
      accrued_allowance <- accrued_allowance + allowance
    }
    
    if(dataset[i,'eom']){
      intrest <- loan_balance * interest_rate_monthly
      accrued_interest <- accrued_interest + intrest
      loan_balance <- loan_balance + intrest
    }
    
    asset_balance <- asset_balance * (1 + dataset[i,'ror'])
    asset_balance_low <- asset_balance * (1 + dataset[i,'ror_low'])
    
    equity_low = asset_balance_low - loan_balance
    
    if(equity_low <= 0){
      loan_balance <- loan_balance - asset_balance_low
      asset_balance <- 0
      dataset[i, 'margin_call'] <- TRUE
    }
    
    dataset[i,'asset_balance'] <- asset_balance
    dataset[i,'loan_balance'] <- loan_balance
    dataset[i,'accrued_interest'] <- accrued_interest
    dataset[i,'accrued_allowance'] <- accrued_allowance
    
  }
  
  dataset$equity <- dataset$asset_balance - dataset$loan_balance
  dataset$cumulative_profit <- dataset$equity - dataset$accrued_allowance
  return(dataset)
}

run_simulations <- function(
    dataset_invest,
    dataset_sample,
    allowance,
    gearing_factor,
    interest_rate,
    p_level,
    sim_count,
    resample,
    simulate){
  library("parallel")
  library("foreach")
  library("doParallel")
  
  core_count <- detectCores() - 2
  cl <- makeCluster(core_count)
  registerDoParallel(cl, cores = core_count)
  
  simres = foreach(i = 1:sim_count, .combine = rbind) %dopar% {
    try({
      dataset_resampled <- resample(dataset_invest, dataset_sample)
      dataset_resampled <- simulate(dataset_resampled, allowance, gearing_factor, interest_rate)
      cumulative_profit <- dataset_resampled$cumulative_profit
      margin_called <- any(dataset_resampled$margin_call)
      
      final_entry <- tail(dataset_resampled,1)
        
      final_cumulative_return <- final_entry$cumulative_profit / final_entry$accrued_allowance
      final_cumulative_profit <- final_entry$cumulative_profit
      
      ror_mean <- mean(dataset_resampled$ror)
      ror_var <- var(dataset_resampled$ror)
      res <- list(margin_called, cumulative_profit, final_cumulative_return, final_cumulative_profit, ror_mean, ror_var)

      return(res)
    })
  }
  stopCluster(cl)
  
  margin_call_prob <- sum(unlist(simres[,1]))/sim_count
  cumulative_profits <- unlist(simres[,2])
  cumulative_return_sd <- sd(unlist(simres[,3]))
  cumulative_profit_sd <- sd(unlist(simres[,4]))
  ror_mean <- mean(unlist(simres[,5]))
  ror_sd <- sqrt(mean(unlist(simres[,6])))
  
  dataset_invest_size <- nrow(dataset_invest)
  
  profit_matrix <- matrix(cumulative_profits, ncol = dataset_invest_size, byrow = TRUE)

  p_lower <- (1-p_level)/2
  p_upper <- 1-p_lower
  p_interval <- c(p_lower, 0.5, p_upper)
  
  bands <- data.frame(t(apply(profit_matrix, 2, quantile, p_interval, na.rm=TRUE)))
  names(bands) <- c('cumulative_profit_lower', 'cumulative_profit_mid', 'cumulative_profit_upper')
  bands$quote_date <- dataset_invest$quote_date
  
  return(list(bands, margin_call_prob, cumulative_return_sd, cumulative_profit_sd, ror_mean, ror_sd))
}

# Run scenarios

for(i in 1:nrow(scenarios)){
  scenario <- scenarios[i,]
  scenario_id <- as.character(scenario$id)
  dataset_sample <- fetch_dataset(scenario$sample_start, scenario$sample_end, exchange, instrument)
  dataset_invest <- fetch_dataset(scenario$invest_start, scenario$invest_end, exchange, instrument)
  
  simres <- run_simulations(
    dataset_invest, 
    dataset_sample, 
    allowance, 
    scenario$gearing_factor, 
    scenario$interest_rate, 
    p_level, 
    sim_count, 
    resample, 
    simulate)
  bands <- simres[[1]]
  margin_call_prob <- simres[[2]]
  cumulative_return_sim_sd <- simres[[3]]
  cumulative_profit_sim_sd <- simres[[4]]
  daily_ror_mean <- simres[[5]]
  daily_ror_sd <- simres[[6]]
  
  dataset <- simulate(dataset_invest, allowance, scenario$gearing_factor, scenario$interest_rate)
  dataset <- merge(x=dataset, y=bands, by='quote_date', all.x=TRUE)
  
  days <- nrow(dataset)
  
  dataset$cumulative_return <- dataset$cumulative_profit / dataset$accrued_allowance
  dataset$cumulative_return_lower <- dataset$cumulative_profit_lower / dataset$accrued_allowance
  dataset$cumulative_return_mid <- dataset$cumulative_profit_mid / dataset$accrued_allowance
  dataset$cumulative_return_upper <- dataset$cumulative_profit_upper / dataset$accrued_allowance
  
  daily_risk_free_return <- (1+risk_free_return) ^ (1/365) - 1 
  #365 days becuase we are including all days of the year in the ror_mean and ror_sd calculations.
  daily_sharp_ratio <- (daily_ror_mean - daily_risk_free_return) / daily_ror_sd
  
  # Output
  
  last_entry <- tail(dataset,1)
  
  print(scenario_id)
  print(paste0('<tr><td>Days</td><td>',days,'</td></tr>'))
  print(paste0('<tr><td>Gearing Factor</td><td>',format(round(scenario$gearing_factor,1),nsmall=1),'x</td></tr>'))
  print(paste0('<tr><td>Interest Rate</td><td>',format(round(scenario$interest_rate*100,2),nsmall=2),'%</td></tr>'))
  print(paste0('<tr><td>Crash</td><td>',substring(scenario_id, nchar(scenario_id)) == "C",'</td></tr>'))
  print(paste0('<tr><td>Estimated Probability of Margin Call</td><td>',format(round(margin_call_prob*100,2),nsmall=2),'%</td></tr>'))
  print(paste0('<tr><td>Cum. Return - Standard Deviation</td><td>',format(round(cumulative_return_sim_sd*100,2),nsmall=2),'%</td></tr>'))
  print(paste0('<tr><td>Cum. Return - Lower (5%)</td><td>',format(round(last_entry$cumulative_return_lower*100,2),nsmall=2),'%</td></tr>'))
  print(paste0('<tr><td>Cum. Return - Mid/Expected (50%)</td><td>',format(round(last_entry$cumulative_return_mid*100,2),nsmall=2),'%</td></tr>'))
  print(paste0('<tr><td>Cum. Return - Upper (95%)</td><td>',format(round(last_entry$cumulative_return_upper*100,2),nsmall=2),'%</td></tr>'))
  print(paste0('<tr><td>Sharp Ratio</td><td>',format(round(daily_sharp_ratio,4),nsmall=4),'</td></tr>'))
  print(paste0('<tr><td>Cum. Profit - Standard Deviation</td><td>NOK ',round(cumulative_profit_sim_sd,0),'</td></tr>'))
  print(paste0('<tr><td>Cum. Profit - Lower (5%)</td><td>NOK ',round(last_entry$cumulative_profit_lower,0),'</td></tr>'))
  print(paste0('<tr><td>Cum. Profit - Mid/Expected (50%)</td><td>NOK ',round(last_entry$cumulative_profit_mid,0),'</td></tr>'))
  print(paste0('<tr><td>Cum. Profit - Upper (95%)</td><td>NOK ',round(last_entry$cumulative_profit_upper,0),'</td></tr>'))
  print(paste0('<tr><td>Accrued Allowance</td><td>NOK ',round(last_entry$accrued_allowance,0),'</td></tr>'))

  library(scales)
  library(ggplot2)
  library(grid)
  library(gridExtra)
  library(gtable)
  
  ret_plot <- ggplot() +
    geom_line(data = dataset, aes(x=quote_date, y=cumulative_return, colour = "Actual"), size=1) +
    geom_line(data = dataset, aes(x=quote_date, y=cumulative_return_lower, colour = "Lower (5%)"), size=1) +
    geom_line(data = dataset, aes(x=quote_date, y=cumulative_return_mid, colour = "Mid (50%)"), size=1) +
    geom_line(data = dataset, aes(x=quote_date, y=cumulative_return_upper, colour = "Upper (95%)"), size=1) +
    labs(x="Time", y="Cumulative Return", colour="") +
    scale_y_continuous(labels = percent, breaks = seq(scenario$ret_low, scenario$ret_high, by = scenario$ret_step), limits = c(scenario$ret_low, scenario$ret_high)) +
    scale_x_date(date_breaks = "2 month", date_labels="%b %y") + 
    theme(legend.position="bottom") +
    ggtitle(paste0(scenario_id,": Cumulative Returns"))
  
  prof_plot <- ggplot() +
    geom_line(data = dataset, aes(x=quote_date, y=cumulative_profit, colour = "Actual"), size=1) +
    geom_line(data = dataset, aes(x=quote_date, y=cumulative_profit_lower, colour = "Lower (5%)"), size=1) +
    geom_line(data = dataset, aes(x=quote_date, y=cumulative_profit_mid, colour = "Mid (50%)"), size=1) +
    geom_line(data = dataset, aes(x=quote_date, y=cumulative_profit_upper, colour = "Upper (95%)"), size=1) +
    labs(x="Time", y="Profit", colour="") +
    scale_y_continuous(labels = comma, breaks = seq(scenario$prof_low, scenario$prof_high, by = scenario$prof_step), limits = c(scenario$prof_low, scenario$prof_high)) +
    scale_x_date(date_breaks = "2 month", date_labels="%b %y") + 
    theme(legend.position="bottom") +
    ggtitle(paste0(scenario_id,": Profit & Loss"))
  
  bal_plot <- ggplot() +
    geom_line(data = dataset, aes(x=quote_date, y=accrued_allowance, colour = "Accrued Allowance"), size=1) +
    geom_line(data = dataset, aes(x=quote_date, y=accrued_allowance+cumulative_profit_lower, colour = "Assets Lower (5%)"), size=1) +
    geom_line(data = dataset, aes(x=quote_date, y=accrued_allowance+cumulative_profit_mid, colour = "Assets Mid (50%)"), size=1) +
    geom_line(data = dataset, aes(x=quote_date, y=accrued_allowance+cumulative_profit_upper, colour = "Assets Upper (95%)"), size=1) +
    labs(x="Time", y="Balance", colour="") +
    scale_y_continuous(labels = comma, breaks = seq(scenario$bal_low, scenario$bal_high, by = scenario$bal_step), limits = c(scenario$bal_low, scenario$bal_high)) +
    scale_x_date(date_breaks = "2 month", date_labels="%b %y") + 
    theme(legend.position="bottom") +
    ggtitle(paste0(scenario_id,": Balances"))
  
  ret_grob <- ggplotGrob(ret_plot)
  prof_grob <- ggplotGrob(prof_plot)
  bal_grob <- ggplotGrob(bal_plot)
  
  maxWidth = grid::unit.pmax(ret_grob$widths[2:5], prof_grob$widths[2:5], bal_grob$widths[2:5])
  
  ret_grob$widths[2:5] <- maxWidth
  prof_grob$widths[2:5] <- maxWidth
  bal_grob$widths[2:5] <- maxWidth
  
  ggsave(filename = paste0('output/return_', i, '_', scenario_id, '.png'), ret_grob)
  ggsave(filename = paste0('output/profit_', i, '_', scenario_id, '.png'), prof_grob)
  ggsave(filename = paste0('output/balance_', i, '_', scenario_id, '.png'), bal_grob)
}