######################################################
# Stat 6721 - Assignment 4
#
# Sarah Rathwell 
#
# Objective - Investigate Merck and Pfizer for and
#             trade on 2015 - Model selection cont.
######################################################
rm(list=ls())
library(quantmod)
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: TTR
## Version 0.4-0 included new data defaults. See ?getSymbols.
######################################################
cat(" Fourth look - Merk vs Pfizer.\n\n",
    "Last modification:",date(),'\n')
##  Fourth look - Merk vs Pfizer.
## 
##  Last modification: Fri Nov 18 13:42:23 2016
######################################################


MERK <- getSymbols('MRK', from = '2010-01-01',to = '2014-12-31', 
                   adjust = T, auto.assign = FALSE)
##     As of 0.4-0, 'getSymbols' uses env=parent.frame() and
##  auto.assign=TRUE by default.
## 
##  This  behavior  will be  phased out in 0.5-0  when the call  will
##  default to use auto.assign=FALSE. getOption("getSymbols.env") and 
##  getOptions("getSymbols.auto.assign") are now checked for alternate defaults
## 
##  This message is shown once per session and may be disabled by setting 
##  options("getSymbols.warning4.0"=FALSE). See ?getSymbols for more details.
PFE <- getSymbols('PFE', from = '2010-01-01',to = '2014-12-31', 
                  adjust = T, auto.assign = FALSE)
head(MERK)
##            MRK.Open MRK.High  MRK.Low MRK.Close MRK.Volume MRK.Adjusted
## 2010-01-04 30.24979 30.60306 30.02797  30.40589   13896500     28.72234
## 2010-01-05 30.64414 30.76738 30.34016  30.52912   14744800     28.83875
## 2010-01-06 30.57020 30.99741 30.31552  30.93990   15239900     29.22679
## 2010-01-07 30.78381 31.13708 30.72630  30.98920   11916600     29.27335
## 2010-01-08 31.11243 31.12886 30.75094  30.97277   10767100     29.25783
## 2010-01-11 31.12065 31.15351 30.80845  31.09600    9581800     29.37424
head(PFE)
##            PFE.Open PFE.High  PFE.Low PFE.Close PFE.Volume PFE.Adjusted
## 2010-01-04 15.13582 15.69088 15.11097  15.68260   52086000     14.62283
## 2010-01-05 15.67431 15.68260 15.36778  15.45892   43372800     14.41426
## 2010-01-06 15.45892 15.58318 15.33465  15.40921   41405100     14.36792
## 2010-01-07 15.44235 15.46720 15.29322  15.35122   39428000     14.31384
## 2010-01-08 15.42578 15.50034 15.34293  15.47548   30407700     14.42971
## 2010-01-11 15.59975 15.69917 15.46720  15.59975   32447500     14.54558
MRK.a <- MERK[,6]
PFE.a <- PFE[,6]

MRK.ld <- diff(log(MRK.a))
MRK.ld <- MRK.ld[-1,]
PFE.ld <- diff(log(PFE.a))
PFE.ld <- PFE.ld[-1,]


PFE.arma10 <- arima(PFE.ld, order=c(1,0,0))
PFE.arma01 <- arima(PFE.ld, order=c(0,0,1))

# --> 1 parameter


MRK.arma50 <- arima(MRK.ld, order=c(5,0,0))
MRK.arma32 <- arima(MRK.ld, order=c(3,0,2))
MRK.arma23 <- arima(MRK.ld, order=c(2,0,3))
MRK.arma05 <- arima(MRK.ld, order=c(0,0,5))

# --> 5 parameters


par(mfrow=c(1,1), oma=c(0,0,2,0))

acf(ts.union(ts(PFE.ld),ts(MRK.ld)))
mtext('Cross Correlation for MRK ~ PFE', outer = TRUE, cex = 1.5)

acf(ts.union(ts(MRK.ld),ts(PFE.ld)))
mtext('Cross Correlation for PFE ~ MRK', outer = TRUE, cex = 1.5)

# run PFE = MRK
# m2 = 0, m1 = 5 (or 1)

union <- ts.union(ts(MRK.ld), ts(PFE.ld))
colnames(union) <- c("MRK", "PFE")
acf(union)
mtext('Cross Correlation for PFE ~ MRK', outer = TRUE, cex = 1.5)

par(mfrow=c(1,2), oma=c(0,0,0,0))
cpgram(PFE.ld, main= "PFE")
cpgram(MRK.ld, main="MRK")

par(mfrow=c(1,2), oma=c(0,0,0,0))
cpgram(PFE.arma10$residuals, main="PFE: arma 10")
cpgram(PFE.arma01$residuals, main="PFE: arma 01")

par(mfrow=c(2,2), oma=c(0,0,0,0))
cpgram(MRK.arma50$residuals, main="MRK: arma 50")
cpgram(MRK.arma32$residuals, main="MRK: arma 32")
cpgram(MRK.arma23$residuals, main="MRK: arma 23")
cpgram(MRK.arma05$residuals, main="MRK: arma 05")

paste(PFE.arma10$aic, PFE.arma01$aic)
## [1] "-7620.41436369809 -7620.36907928162"
# AR1 is marginally smaller

paste(MRK.arma50$aic, MRK.arma32$aic, MRK.arma23$aic, MRK.arma05$aic)
## [1] "-7566.06326450209 -7561.14593501032 -7561.15837102485 -7564.00913126507"
# AR5 is marginally smaller

par(mfrow=c(1,2), oma=c(0,0,0,0))
cpgram(PFE.arma10$residuals, main="PFE: AR(1)")
cpgram(MRK.arma50$residuals, main="MRK: AR(5)")