Create Arima Code

# Please Install these two Libraries
library(Quandl);
library(xlsx);
library(tseries, quietly = T)
library(timeSeries)
library(forecast)
library(TTR)
library(quantmod)
library(fArma)
library(MASS)
library(forecast, quietly = T)
library(outliers)
require(quantmod);
library(lattice)
library(PerformanceAnalytics)
library(TSA)
library(gridExtra)

#Get the data
Symbol = "SPY"
getSymbols(Symbol,  src = "yahoo",from = as.Date("2009-01-10"),to = as.Date("2013-07-31"));
RawData<- SPY
StockData = RawData

#Get the Close data
close1<-(StockData$SPY.Close)
#Plot the time series of the close data
plot(close1, ylab = "Closing Price", main = "Closing Price Time Series")
lines(close1, ylab = "Closing Price",col='red', main = "Closing Price Time Series")

#Plot the Autocorrelation Plot for the close, check stationarity and order of time series
acfPlot(close1)
kpss.test(close1)
auto.arima(close1)

#differntiating the close time series
ret2 <- 100*diff(log(close1))
ret2 [is.na(ret2 )] <- 0 

#stationarity test and finding the order of the time series
kpss.test(ret2)
auto.arima(ret2)

#Plot the Time series plot of the differentiated Close Series (Returns) [Box-COx if needed]
plot(ret2, ylab = "Log Daily Returns"  ,main = "Log Daily Returns Time Series")
lines(ret2, ylab = "Log Daily Returns" ,col='purple' ,main = "Log Daily Returns Time Series")
windows()
#BoxCox.lambda(ret2, method=c("guerrero","loglik"), lower=0, upper=2)
#ret2 <- BoxCox(ret2 ,BoxCox.lambda(ret2, method=c("guerrero","loglik"), lower=0, upper=2))

#Time series Plot and autocorrelation plot of the returns
plot(ret2, ylab = "Log Daily Returns",  main = "Log Daily Returns Time Series")
lines(ret2, ylab = "Log Daily Returns" ,col='purple' ,main = "Log Daily Returns Time Series")
windows()
ret2[is.na(ret2)] <- 0
acfPlot(ret2)
ret1<-ret2

#Removing the outliers if any and finding order and plotting the time series
rm.outlier(ret1, fill = FALSE, median = FALSE, opposite = FALSE)
windows()
plot(ret1, ylab = "Log Daily Returns"  ,main = "Log Daily Returns Time Series")
lines(ret1, ylab = "Log Daily Returns",col='purple',  main = "Log Daily Returns Time Series")
kpss.test(ret1)
auto.arima(ret1)

#Converting the time series into vector for analysis
ret1<-as.numeric(ret1)
auto.arima(ret1)


Length7<-length(ret1)-50
Length8<-Length7/length(ret1)
ret_train3 <-ret1[1:(Length8* length(ret1))]  # Train dataset
ret_test3 <- ret1[((Length8* length(ret1)) + 1)]  # Test dataset
auto.arima(ret_train3)

#Cross Validation - Train, Test and Predict, and creating data frame with sign indicators and trend indicators and checking accuracy of predicion
j=1
#creating data frame for actual(+1 and -1) and for predicted ones and also backtesting on the 50 days window (Cross Validation)
sign2=NULL;
sign3=NULL;
mult<-vector()
arma.preds4<-vector()
ret_test4<-vector()
for(i in 50:1){
  Length7<-length(ret1)-i
  Length8<-Length7/length(ret1)
  
  ret_train3 <-ret1[1:(Length8* length(ret1))]  # Train dataset
  ret_test3 <- ret1[((Length8* length(ret1)) + 1)]  # Test dataset
  
  #Predict one day ahead over 50 days from last date of the data
  fit3 <- arima(ret_train3, order = c(1, 0, 1))
  arma.preds3 <- predict(fit3, n.ahead = 1,type="prob")$pred
  
  
  arma.preds4[j]<-arma.preds3
  ret_test4[j]<-ret_test3
  
 #FInding out the Trend indicator
  if(arma.preds3>ret_train3[length(ret_train3)]){
    sign2[j]=1;
  }
  else{
    sign2[j]=-1;
  }
  
  if(ret_test3>ret_train3[length(ret_train3)]){
    sign3[j]=1;
  }
  else{
    sign3[j]=-1;
  }
  
  mult[j]<-sign2[j]*sign3[j]
  j=j+1
}

#Data Frame Creation with Trend indicator
Testframe<-data.frame(arma.preds4,ret_test4,sign2,sign3,mult)
length(mult[mult==1])

#New Test Frame1 with Sign Indicator
Indicator=vector()
for(i in 1:50)
{
  if(arma.preds4[i]<0)
  {
    Indicator[i]=-1
  }
  else{
    Indicator[i]=1
  }
}
retfortest<-ret1[(length(ret1)-49):length(ret1)]/100
#New Test Frame
Testframe1<-data.frame(arma.preds4,sign2,retfortest,sign3,mult, Indicator)


#Testing the Arma movement
sign2<-Testframe1$sign2;
Retfortest<-Testframe1$retfortest;
Armagrowth=vector()
Armagrowth1=vector()
Armagrowth1[1]=1000

#Arima vs Buy Hold comparison with our performance factors as the keys (Trend and sign Indicator)
for(i in 1:50)
{
  if((Indicator[i]==-1|sign2[i]==-1)){
    
    Armagrowth[i]=Armagrowth1[i];
    Armagrowth1[i+1]=Armagrowth[i];
  }
  else
  {
    
    Armagrowth[i]=Armagrowth1[i]+Retfortest[i]*Armagrowth1[i];
    Armagrowth1[i+1]=Armagrowth[i];
  }
}



#Testing the Actual movement
Actualgrowth=vector()
Actualgrowth1=vector()
Actualgrowth1[1]=1000
for(i in 1:50)
{
  Actualgrowth[i]<-Actualgrowth1[i]+Retfortest[i]*Actualgrowth1[i];
  Actualgrowth1[i+1]=Actualgrowth[i];
}

#Time Series of the predicted and actual portfolio value
Armagrowth<- xts(Armagrowth,order.by=index(ret2[(length(ret2)-49):length(ret2)]))
Actualgrowth<- xts(Actualgrowth,order.by=index(ret2[(length(ret2)-49):length(ret2)]))

#Final Plot of Arima annd actual portfolio value from buy and hold strategy
ts.plot(ts(Armagrowth),ts(Actualgrowth),col=c("blue","red"),ylab="Portfolio Value",xlab="Period(2013/05/23-2013/07/31)",main="Comparison of ARIMA vs Buy/Hold for S&P 500 over last day of 50 days")
legend(-0.4,1022,legend=c("ARIMA", "Buy/Hold"), col=c("blue","red"),lty =1)
windows()


#Final Plot of Arima annd actual portfolio value from buy and hold strategy
yrngMin = abs(min(Armagrowth,Actualgrowth))
yrngMax = abs(max(Armagrowth,Actualgrowth))
chartSeries(Armagrowth,Actualgrowth, log.scale = TRUE, name="Portfolio Value", yrange=c(yrngMin, yrngMax))
addTA(Actualgrowth,on=1,col="gold")


#Transfer Function Model and cross correlation
CloseT<-(StockData$SPY.Adjusted)
OpenT<-(StockData$SPY.Open)
LowT<-(StockData$SPY.Low)
HighT<-(StockData$SPY.High)
VolumeT<-(StockData$SPY.Volume)


ccf(as.numeric(CloseT),as.numeric(OpenT),ylab="Close and Opening Price Correlation",xlab="Time")
windows()
ccf(as.numeric(CloseT),as.numeric(LowT),ylab="Close and Low Price Correlation",xlab="Time")
windows()
ccf(as.numeric(CloseT),as.numeric(HighT),ylab="Close and High Price Correlation",xlab="Time")
windows()
ccf(as.numeric(CloseT),as.numeric(VolumeT),ylab="Close and Volume Price Correlation",xlab="Time")

#kpss test for all
kpss.test(CloseT)
kpss.test(OpenT)
kpss.test(LowT)
kpss.test(HighT)
kpss.test(VolumeT)

#Time series plotting
ts.plot(CloseT)
windows()
ts.plot(VolumeT)
windows()

#Process of making the time series stationary and checking order [Plotting the time series]
ret4 <- 100*diff(log(CloseT))
ret4[is.na(ret4)] <- 0 
ret5<-as.numeric(ret4)
plot(ret5)
windows()
auto.arima(ret5)

volumediff <- 100*diff(log(VolumeT))
volumediff[is.na(volumediff)] <- 0 
rm.outlier(volumediff, fill = FALSE, median = FALSE, opposite = FALSE)
kpss.test(volumediff)
plot(volumediff)
windows()
volumediff<-as.numeric(volumediff)
auto.arima(volumediff)

#ret of open
openret <- 100*diff(log(OpenT))
openret[is.na(openret)] <- 0 
plot(openret)
windows()
openret<-as.numeric(openret)
auto.arima(openret)


#ret of high
highret <- 100*diff(log(HighT))
highret[is.na(highret)] <- 0 
plot(highret)
windows()
highret<-as.numeric(highret)
auto.arima(highret)


#ret of low
lowret <- 100*diff(log(LowT))
lowret[is.na(lowret)] <- 0 
ts.plot(lowret)
windows()
lowret<-as.numeric(lowret)
auto.arima(lowret)

#transfer function for prediction
#regParams = ts.union(ts(lowret), ts(highret), ts(openret))
#model10 = arimax(ret5,order=c(1,0,1),xreg=regParams)


#selecting the proper model for  transfer function using ridge regression, lasso, stepwise and elastic nets
Selectframe<-data.frame(ret1,openret,lowret,highret,volumediff)
fitridge<-select(lm.ridge(ret1~(openret+lowret+highret+volumediff), Selectframe))
fitridge1<-(lm.ridge(ret1~(openret+lowret+highret+volumediff), Selectframe))

fitstep1 <- lm(ret1~(openret+lowret+highret+volumediff),data=Selectframe)
stepAIC(fitstep1,direction="both",steps=1000)


#converting the data into vector for our analysis
lowret<-as.numeric(lowret)
highret<-as.numeric(highret)
openret<-as.numeric(openret)
volumediff<-as.numeric(volumediff)
nreg1=vector()
j=1
#creating data frame for actual(+1 and -1) and for predicted ones for transfer function model, cross validation
sign4=NULL;
sign5=NULL;
mult1<-vector()
arma.preds5<-vector()
arma.preds6<-vector()
ret_test5<-vector()
arma.preds6<-vector()
auto.arima
#Transfer Function Implementation for backtesting 50 days length
for(i in 50:1){
  Length7<-length(ret1)-i
  Length8<-Length7/length(ret1)
  
  ret_train3 <-ret1[1:(Length8* length(ret1))]  # Train dataset
  ret_test3 <- ret1[((Length8* length(ret1)) + 1)]  # Test dataset
  
  
  regParams <- as.numeric((lowret))
  #Transfer Function Implementation after selection of the variable and the transfer function order
  xreg1<- regParams[1:(Length8* length(ret1))]
  model10 = arima(ret_train3,order=c(1,0,1),xreg= xreg1);
  #predicting one day ahead over 50 days
  arma.preds5[51-i] <- predict( model10,n.ahead = 1,newxreg= xreg1[((Length8* length(ret1))-4)])
}
arma.preds5<-as.numeric(arma.preds5)

#Finding out the Trend Indicator and Sign Indicator over the required length
ret_test3<-vector()
ret_train3<-vector()
for(i in 50:1){
  j=51-i
  Length7<-length(ret1)-i
  Length8<-Length7/length(ret1)
  ret_train3[51-i] <-ret1[(Length8* length(ret1))]
  ret_test3[51-i] <- ret1[((Length8* length(ret1)) + 1)]
}

sign4=NULL;
sign5=NULL;
mult1<-vector()
for(i in 1:50){
  if(arma.preds5[i]>ret_train3[i]){
    sign4[i]=1;
  }
  else{
    sign4[i]=-1;
  }
  
  if(ret_test3[i]>ret_train3[i]){
    sign5[i]=1;
  }
  else{
    sign5[i]=-1;
  }
  
  mult1[i]<-sign4[i]*sign5[i]
  
  
}

#Data Frame creation of the predicted and actual value (test) with the indicators (sign and trend indicators)
Testframe3<-data.frame(arma.preds5,ret_test3,sign4,sign5,mult1)
length(mult1[mult1==1])

#New Test Frame1
Indicator1=vector()
for(i in 1:50)
{
  if(arma.preds5[i]<0)
  {
    Indicator1[i]=-1
  }
  else{
    Indicator1[i]=1
  }
}
retfortest1<-ret1[(length(ret1)-49):length(ret1)]/100
Testframe4<-data.frame(arma.preds5,sign4,Indicator1,retfortest1)


#Testing the Arma movement (with our performance parameters)

Retfortest<-Testframe1$retfortest;
Armagrowth3=vector()
Armagrowth1=vector()
Armagrowth1[1]=1000


for(i in 1:50)
{
  if((Indicator1[i]==-1|sign4[i]==-1))
  {
    
    
    Armagrowth3[i]=Armagrowth1[i];
    Armagrowth1[i+1]=Armagrowth3[i];
  }
  else
  {
    Armagrowth3[i]=Armagrowth1[i]+Retfortest[i]*Armagrowth1[i];
    Armagrowth1[i+1]=Armagrowth3[i];
  }
}


ts.plot( ts(Armagrowth3))

#Testing the Actual movement
Actualgrowth2=vector()
Actualgrowth3=vector()
Actualgrowth3[1]=1000
for(i in 1:50)
{
  Actualgrowth2[i]<-Actualgrowth3[i]+Retfortest[i]*Actualgrowth3[i];
  Actualgrowth3[i+1]=Actualgrowth2[i];
}
windows()


#Time Series of the predicted and actual portfolio value (Buy and Hold)
Armagrowth3<- xts(Armagrowth3,order.by=index(ret2[(length(ret2)-49):length(ret2)]))
Actualgrowth2<- xts(Actualgrowth2,order.by=index(close1[(length(ret2)-49):length(ret2)]))

windows()
#Final Plot of Arima annd actual portfolio value (Buy and Hold)
ts.plot(ts(Armagrowth3),ts(Actualgrowth2),col=c("blue","red"),ylab="Portfolio Value",xlab="Period(2013/05/23-2013/07/31)",main="Comparison of Transfer Function vs Buy/Hold for S&P 500 over last day of 50 days")
legend(-0.6,1022,legend=c("Transfer Function", "Buy/Hold"), col=c("blue","red"),lty =1)

#Final Plot of Arima actual portfolio value (Buy and Hold)
yrngMin = abs(min(Armagrowth3,Actualgrowth2))
yrngMax = abs(max(Armagrowth3,Actualgrowth2))
chartSeries(Armagrowth3,log.scale = TRUE, name="Portfolio Value", yrange=c(yrngMin, yrngMax))
addTA(Actualgrowth2,on=1,col="gold")


Armagrowth3<-as.vector(Armagrowth3)
Armagrowth<-as.vector(Armagrowth)
Actualgrowth2<-as.vector(Actualgrowth2)
Actualgrowth<-as.vector(Actualgrowth)
#Statistics
#statistics of cumulative returns and sharp ratio for ARIMA
Arimacumulativereturn<-(100*Delt((Armagrowth[1:50])))
Arimacumulativereturn1<-sum(Arimacumulativereturn[2:50])
print(Arimacumulativereturn1)


Actualcumulativereturn<-100*Delt((Actualgrowth[1:50]))
Actualcumulativereturn<-as.numeric(Actualcumulativereturn[2:50])
print(Actualcumulativereturn1<-sum(Actualcumulativereturn))


#statistics of cumulative returns and sharp ratio for Transfer FUnction
Transcumulativereturn<-(100*Delt((Armagrowth3[1:50])))
print(Transcumulativereturn1<-sum(Transcumulativereturn[2:50]))

Trans1cumulativereturn<-100*Delt((Actualgrowth2[1:50]))
print(Trans1cumulativereturn1<-sum(Trans1cumulativereturn[2:50]))

#Sharpe Ratio for Arima
print(Sharperatioactual<-sharpe(Actualgrowth,scale=1))
print(Sharperatioarima<-sharpe(Armagrowth,scale=1))

#Sharpe Ratio for TF
print(Sharperatioactual1<-sharpe(Actualgrowth2,scale=1))
print(Sharperatiotrans<-sharpe(Armagrowth3,scale=1))

#Result Statistics ARIMA
PredAccuracy<-(length(mult[mult==1]))/50
BenchMark = c(Actualcumulativereturn1, Sharperatioactual,PredAccuracy)
Strategy = c(Arimacumulativereturn1, Sharperatioarima,PredAccuracy)
nms = c("Cumulative Return", "Sharpe Ratio","Prediction Accuracy")
result = data.frame(Strategy, BenchMark, row.names = nms)
windows()

accuracytrans<-accuracy(arma.preds5,retfortest1)
accuracyarima<-accuracy(arma.preds4,retfortest)
#Tabulating the Statisics
grid.table((result),gpar.coretext = gpar(fontsize=30,col="slateblue"),gpar.corefill = gpar(col="slateblue",alpha=0.5),h.even.alpha = 0.5)
windows()
grid.table((accuracyarima),gpar.coretext = gpar(fontsize=15,col="slateblue"),gpar.corefill = gpar(col="slateblue",alpha=0.5),h.even.alpha = 0.5)
windows()

#Result Statistics TF
PredAccuracy1<-(length(mult1[mult1==1]))/50
BenchMark = c(Trans1cumulativereturn1, Sharperatioactual1,PredAccuracy1)
Strategy = c(Transcumulativereturn1, Sharperatiotrans,PredAccuracy1)
nms = c("Cumulative Return", "Sharpe Ratio","Prediction Accuracy")
result1 = data.frame(Strategy, BenchMark, row.names = nms)

#Tabulating the Statisics for TF
grid.table((result1),gpar.coretext = gpar(fontsize=30,col="slateblue"),gpar.corefill = gpar(col="slateblue",alpha=0.5),h.even.alpha = 0.5)
windows()
grid.table((accuracytrans),gpar.coretext = gpar(fontsize=15,col="slateblue"),gpar.corefill = gpar(col="slateblue",alpha=0.5),h.even.alpha = 0.5)

#p value tests for ARIMA MOdel
summary(arma(ret1, order = c(1, 1)))