#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

save_plots = FALSE

set.seed(0) 

#
# Data Analysis
#

# Problem 1: EPage 19
#
dat = read.csv("../../BookCode/Data/Stock_FX_Bond.csv",header=TRUE)
names(dat)
attach(dat)
par(mfrow=c(1,2))
plot(GM_AC, type='l')
plot(F_AC, type='l')
par(mfrow=c(1,1))

n = dim(dat)[1]
GMReturn = GM_AC[2:n]/GM_AC[1:(n-1)] - 1
FReturn = F_AC[2:n]/F_AC[1:(n-1)] - 1
par(mfrow=c(1,1))

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter2/chap_2_rlab_prob_1_GM_N_F_returns.eps", onefile=FALSE, horizontal=FALSE) }
plot(GMReturn,FReturn)
grid()
if( save_plots ){ dev.off() }


# Problem 2:
# 
GMLogReturn = log( GM_AC[2:n]/GM_AC[1:(n-1)] )
plot( GMReturn, GMLogReturn )
abline(a=0,b=1)
print( cor( GMReturn, GMLogReturn ) )
print( cor( GMReturn, FReturn ) )


# Problem 3:
#
MSFT_Return = MSFT_AC[2:n]/MSFT_AC[1:(n-1)] - 1
MRK_Return = MRK_AC[2:n]/MRK_AC[1:(n-1)] - 1

print( cor( MSFT_Return, MRK_Return ) )

if( save_plots ){ postscript('../../WriteUp/Graphics/Chapter2/chap_2_rlab_prob_3_MSFT_N_MRK_returns.eps', onefile=FALSE, horizontal=FALSE) }
plot(MSFT_Return, MRK_Return, xlab='MSFT return', ylab='MRK return')
grid()
if( save_plots ){ dev.off() }


# Problem 4: An example MC simultation: 
#
niter = 1.e5 # number of simulations of our price process
below = rep(0,niter) # indice as to whether we have to liquidate our portfolio during this 45 days
set.seed(2009)
for( i in 1:niter ){
  r = rnorm( 45, mean=0.05/253, sd=0.23/sqrt(253) )
  logPrice = log(1.e6) + cumsum(r) # the time series of portfolio values
  minlogP = min(logPrice)
  below[i] = as.numeric(minlogP < log(950000))
}
mean(below)


# Problems 5-7:
#
niter = 1.e5 # number of simulations of our price process
n_days_simulation = 100

profit_target = 1.e6 + 100000.
log_profit_target = log(profit_target) 

stop_loss = 1.e6 - 50000.
log_stop_loss = log(stop_loss) 

above = rep(0,niter) # indices as to whether we have to liquidate our portfolio during this 100 days due to a profit 
below = rep(0,niter) # indices as to whether we have to liquidate our portfolio during this 100 days due to a loss
middle = rep(0,niter) # we don't liquidate our portfolio but hold till the end of the time frame 
profit = rep(0,niter)
strat_return = rep(0,niter) # over the 100 days 
set.seed(2009)
for( i in 1:niter ){
  r = rnorm( n_days_simulation, mean=0.05/253, sd=0.23/sqrt(253) ) # trade for n_days_simulation days 
  logPrice = log(1.e6) + cumsum(r) # the time series of portfolio values
  minlogP = min(logPrice)
  maxlogP = max(logPrice)

  # Find out whether either of the stoploss/profit targets are triggered during this MC price history: 
  #
  sl_triggered = minlogP <= log_stop_loss
  pt_triggered = maxlogP >= log_profit_target

  # only stoploss triggered in our MC simulation of price: 
  # 
  if( sl_triggered & !pt_triggered ){ 
    below[i] = 1
    profit[i] = -50000 # at the moment the value crosses below 950,000 we must exit with a loss of 50K
    days_in_trade = which( logPrice <= log_stop_loss )[1] 
  }

  # only profit target triggered in our MC simulation of price: 
  # 
  if( !sl_triggered & pt_triggered ){ 
    above[i] = 1
    profit[i] = +100000 # at the moment the value crosses above 1,100,000 we exit with a profit of 100K
    days_in_trade = which( logPrice >= log_profit_target )[1] 
  }

  # both values triggered ... need to find out which one happend first
  #
  if( sl_triggered & pt_triggered ){ 
    days_in_trade_to_min = which( logPrice <= log_stop_loss )[1] 
    days_in_trade_to_max = which( logPrice >= log_profit_target )[1] 

    if( days_in_trade_to_min < days_in_trade_to_max ){ # Minimum happens first and we hit our stoploss: 
      below[i] = 1
      profit[i] = -50000
      days_in_trade = days_in_trade_to_min
    }else{                                             # Maximum happens first and we hit our profit target: 
      above[i] = 1 
      profit[i] = +100000
      days_in_trade = days_in_trade_to_max
    }
  }

  # neither values triggered: 
  if( !sl_triggered & !pt_triggered ){
    middle[i] = 1
    profit[i] = exp(logPrice[n_days_simulation]) - 1.e6 # we get how much we are above / below 1
    days_in_trade = n_days_simulation
  }

  # Compute returns: 
  strat_return[i] = profit[i] / 50000.
  strat_return[i] = strat_return[i] / days_in_trade # convert the raw return into daily returns 
}
stopifnot( abs( mean(above)+mean(below)+mean(middle) - 1) < 1.e-6 )
print(mean(above))
print(mean(below))
print(mean(profit))
print(mean(strat_return))


# Problems 8 - 11:
#
set.seed(2012)
n = 253
par(mfrow=c(3, 3))
for( i in (1:9) ){
    logr = rnorm(n, 0.05 / 253, 0.2 / sqrt(253))
    price = c(120, 120*exp(cumsum(logr)))
    plot(price, type='b')
}
par(mfrow=c(1, 1))


# Problems 12 - 13:
#
data = read.csv('../../BookCode/Data/MCD_PriceDaily.csv')
head(data)
adjPrice = data[, 7]

n = length(adjPrice)
rets = adjPrice[2:n] / adjPrice[1:(n-1)] - 1
lrets = log( adjPrice[2:n] / adjPrice[1:(n-1)] )
plot(rets, lrets, type='p', xlab='return', ylab='log return')
grid()
print(cor(rets, lrets))

print(sprintf('returns: mean= %f; std= %f', mean(rets), sd(rets)))
print(sprintf('log-returns: mean= %f; std= %f', mean(lrets), sd(lrets)))

t.test(x=rets, y=lrets, paired=TRUE)


# Problem 16:
#
set.seed(2015)

mu = mean(lrets)
s = sd(lrets)

p = 93.07
T = 20
p_limit = 85.00
p_limit_2 = 84.50
n_mc = 10000
profit_game_1 = rep( NA, n_mc )
profit_game_2 = rep( NA, n_mc )
for( mci in 1:n_mc ){
    logr = rnorm(T, mu, s)
    prices = c(p, p*exp(cumsum(logr)))

    # Play the first game:
    #
    if( min(prices) < p_limit ){
        profit_game_1[mci] = +100
    }else{
        profit_game_1[mci] = -1
    }

    # Play the second game:
    #
    if( min(prices) < p_limit_2 ){
        profit_game_2[mci] = +100+25
    }else{
    if( min(prices) < p_limit ){
        profit_game_2[mci] = +100
    }else{
        profit_game_2[mci] = -1
    }
    }
}
print(sprintf('mean(profit_game_1)= %f', mean(profit_game_1)))
print(sprintf('mean(profit_game_2)= %f', mean(profit_game_2)))