#
# 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 = F

library(urca)

set.seed(0) 

#
# Cointegration of Midcap Prices 
#

# P1 EPage 420
# 
midcapD.ts = read.csv('../../BookCode/Data/midcapD.csv', header=TRUE)
x = midcapD.ts[,2:11] 
prices = exp(apply(x,2,cumsum))
options(digits=3)
summary(ca.jo(prices))

#
# Cointegration Analysis of Yields:
# P2-3 EPage 421
# 
mk.zero2 = read.csv('../../BookCode/Data/mk.zero2.csv', header=TRUE)
mk.maturity = read.csv('../../BookCode/Data/mk.maturity.csv', header=TRUE)
print(mk.maturity$Maturity[2:11] * 12)
summary(ca.jo(mk.zero2[,2:11]))


#
# Cointegration Analysis of Daily Stock Prices:
# P 4 Epage 461
#
CokePepsi = read.table('../../BookCode/Data/CokePepsi.csv', header=TRUE)
ts.plot(CokePepsi)

ts.plot(CokePepsi[, 2] - CokePepsi[, 1])

summary(ca.jo(CokePepsi))


Stock_FX_Bond = read.csv('../../BookCode/Data/Stock_FX_Bond.csv', header=TRUE)
adjClose = Stock_FX_Bond[, seq(from=3, to=21, by=2) ]
ts.plot(adjClose)
summary(ca.jo(adjClose))

summary(ca.jo(adjClose, K=8))

#
# Simulation 
#

# P4-7 EPage 421
# 
mu = 1030000
phi = 0.99

niter = 1.e4 # number of simulations of our price process

profit_target = 1.e6 + 20000.

stop_loss = 1.e6 - 50000.

above = rep(0,niter) # indices as to whether we have to liquidate our portfolio during our trading period 
below = rep(0,niter) # indices as to whether we have to liquidate our portfolio during our trading period 
profit = rep(0,niter) # dollar profit made / lost on this run 
strat_holding_period = rep(0,niter)
set.seed(2009)
for( i in 1:niter ){
  t = 0
  P_t = 1.e6 # initial value of 1M

  # Walk forward in time until we liquidate 
  while( TRUE ){
    # Simulate the portfolios new value: 
    P_tp1 = mu + phi * ( P_t - mu ) + rnorm( 1, mean=0, sd=5000 )
    t = t+1 

    if( P_tp1 >= profit_target ){
      above[i] = 1
      profit[i] = 20000.
      strat_holding_period[i] = t
      break
    }
    if( P_tp1 <= stop_loss ){
      below[i] = 1
      profit[i] = -50000.
      strat_holding_period[i] = t
      break
    }
    
    P_t = P_tp1 
  }
}
strat_return = ( profit / 50000. ) / strat_holding_period # compute daily return
stopifnot( abs( mean(above)+mean(below) - 1) < 1.e-6 )

print(mean(profit))
print(mean(below))
print(mean(strat_holding_period))
print(mean(strat_return) * 253) # expected yearly return