MST

library(quantmod)
library(rugarch)
library("PerformanceAnalytics")
 
getSymbols("SPY", from="1900-01-01")
getSymbols("^GDAXI", from="1900-01-01")
spyRets = na.omit( ROC( Cl( GDAXI) ) )

getSymbols("^GSPC", from="1900-01-01")
spyRets = na.omit( ROC( Cl( GSPC) ) )

getSymbols("IBM", from="1900-01-01")
spyRets = na.omit( ROC( Cl( IBM) ) )

plot(cumsum(spyRets["2008"]))
plot(cumsum(spyRets["2014"]))

spyRets = na.trim( ROC( Cl( SPY ) ) )
#spyRets = na.trim( diff( Cl( SPY ) ) )

# Train over 2000-2004, forecast 2005
ss = spyRets["2006/2008"]
outOfSample = NROW(ss["2008"])

ss = spyRets["2011/2014"]
outOfSample = NROW(ss["2014"])

spec = ugarchspec(
            variance.model=list(garchOrder=c(1,1)),
            mean.model=list(armaOrder=c(1,1), include.mean=T),
            distribution.model="sged")

fit = ugarchfit(spec=spec, data=ss, out.sample=outOfSample)
fore = ugarchforecast(fit, n.ahead=1, n.roll=outOfSample)

# Build some sort of indicator base on the forecasts
ind = xts(head(as.array(fore)[,2,],-1), order.by=index(ss["2008"]))
ind = xts(head(as.array(fore)[,2,],-1), order.by=index(ss["2014"]))

ind = ifelse(ind < 0, 1, -1)
 
# Compute the performance
mm = merge( ss["2008"], ind, all=F )
mod <- mm[,1]*mm[,2]
charts.PerformanceSummary(merge(spyRets["2008"],mod))

mm = merge( ss["2014"], ind, all=F )
mod <- mm[,1]*mm[,2]
charts.PerformanceSummary(merge(spyRets["2014"],mod))

plot(cumsum(mod))

cumprod(1+1*2)

tail(cumprod(mm[,1]*mm[,2]+1))
plot(cumprod(mm[,1]*mm[,2]+1))
# Output (last line): 2005-12-30  1.129232


library(quantmod)
library(fArma)

# Get S&P 500
getSymbols( "^GSPC", from="2000-01-01" )

# Compute the daily returns
GSPC.rets = diff(log(Cl(GSPC)))

# Use only the last two years of returns
GSPC.tail = as.ts( tail( GSPC.rets, 500 ) )

# Fit the model
GSPC.arma = armaFit( formula=~arma(2,2), data=GSPC.tail )


xxArma = armaFit( xx ~ arma( 5, 1 ), data=xx )
xxArma@fit$aic

findBestArma = function( xx, minOrder=c(0,0), maxOrder=c(5,5), trace=FALSE )
{
   bestAic = 1e9 
   len = NROW( xx )
   for( p in minOrder[1]:maxOrder[1] ) for( q in minOrder[2]:maxOrder[2] )
   {   
      if( p == 0 && q == 0 ) 
      {   
         next
      }   
 
      formula = as.formula( paste( sep="", "xx ~ arma(", p, ",", q, ")" ) ) 
 
      fit = tryCatch( armaFit( formula, data=xx ),
                      error=function( err ) FALSE,
                      warning=function( warn ) FALSE )
      if( !is.logical( fit ) ) 
      {   
         fitAic = fit@fit$aic
         if( fitAic < bestAic )
         {   
            bestAic = fitAic
            bestFit = fit 
            bestModel = c( p, q ) 
         }   
 
         if( trace )
         {   
            ss = paste( sep="", "(", p, ",", q, "): AIC = ", fitAic )
            print( ss )
         }   
      }   
      else
      {   
         if( trace )
         {   
            ss = paste( sep="", "(", p, ",", q, "): None" )
            print( ss )
         }   
      }   
   }   
 
   if( bestAic < 1e9 )
   {   
      return( list( aic=bestAic, fit=bestFit, model=bestModel ) ) 
   }   
 
   return( FALSE )
}

library(quantmod)
library(fArma)
 
getSymbols("SPY", from="1900-01-01")
SPY.rets = diff(log(Ad(SPY)))
SPY.arma = armaFit(~arma(0, 2), data=as.ts(tail(SPY.rets,500)))
predict(SPY.arma, n.ahead=1, doplot=F)

# Now, to build an indicator for back testing, one can walk the 
# daily return series and at each point perform the steps we covered so 
# far. The main loop looks like (in pseudocode):

for(ii in history:length(dailyRetSeries))
{
   tt = as.ts(tail(head(dailyRetSeries, ii), history))
   ttArma = findBestArma()
   predict(ttArma, n.ahead=1, doplot=F)
}

library(quantmod)
library(fGarch)
 
getSymbols("SPY", from="1900-01-01")
SPY.rets = diff(log(Ad(SPY)))
SPY.garch = garchFit(~arma(0, 2) + garch(1, 1), data=as.ts(tail(SPY.rets, 500)))
predict(SPY.garch, n.ahead=1, doplot=F)