## MA models
 
# Signal to noise example:

epsilon=rnorm(500,sd=1)
eta=rnorm(500,sd=sqrt(0.5))
z=cumsum(eta)
y=z+epsilon 

Delta.y=diff(y)

plot(y,type="l",main="Signal plus noise",ylab="y")
plot(Delta.y,type="l",main=expression(paste(1^{st}," differences")),ylab=expression(Delta * y))
## for more details about the labels type ?plotmath


m1=arima(Delta.y,order=c(0,0,1))
m1        
> m1

Call:
arima(x = Delta.y, order = c(0, 0, 1))

Coefficients:
          ma1  intercept
      -0.5265     0.0148             ## The model is Delta.y(t) =  0.0148+a(t)-0.5265a(t-1).
s.e.   0.0391     0.0290


tsdiag(m1)
predict(m1,5)

acf(Delta.y,lag.max=10)
pacf(Delta.y,lag=10)
acf(Delta.y,lag.max=10)$acf


## Example overlapping / nonoverlapping returns.
sp500monthly00.09=read.csv("GSPCmonthly2000-May2009.csv",header=T)
sp500monthly00.09.price= sp500monthly00.09$Adj.Close[length(sp500monthly00.09$Adj.Close):1] ## 113 months

library(fBasics)
#TO calculate returns automatically:
returns(sp500monthly00.09.price)  ## default gives continuously compunded returns. Note the first obs is NA
returns(sp500monthly00.09.price,method="compound")  ## same as above

returns( sp500monthly00.09.price,method="simple")  ## gives the simple returns  
returns( sp500monthly00.09.price,method="discrete")  ## same as simple returns above

sp500monthly00.09.return=  returns( sp500monthly00.09.price)[-1] ## to eliminate the first NA element
## To calculate overlapping returns look at documentation:
?RollingAnalysis

library(fTrading)
# we need to apply sum function:
sp500overlapping.returns= rollFun(sp500monthly00.09.return,n=12,FUN=sum)

acf(sp500overlapping.returns)  ## Note the MA(10) behavior
pacf(sp500overlapping.returns) ## Interesting the spike at lag 13

eacf.tables(sp500overlapping.returns) # see code for eacf

## Another way: find the lowest AIC
max.p=15
max.q=15
model.aic=matrix(,nrow=max.p+1,ncol=max.q+1)
for (p in 0:max.p){
  for (q in 0:max.q){
  model.aic[p+1,q+1]= arima(sp500overlapping.returns,order=c(p,0,q))$aic
  }
}                   ## takes forever

Error in arima(sp500overlapping.returns, order = c(p, 0, q)) : 
  non-stationary AR part from CSS
In addition: There were 15 warnings (use warnings() to see them)

print(model.aic)

            [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]      [,8]      [,9]     [,10]     [,11]     [,12]     [,13]     [,14]     [,15]     [,16]
 [1,]  -33.48219 -145.4638 -189.3015 -223.1889 -237.5765 -254.3533 -258.3738 -261.5222 -261.7785 -267.5060 -266.2108 -289.7048 -291.7007 -289.7155 -288.1641 -286.4500
 [2,] -275.02779 -275.3925 -273.7859 -273.9294 -271.9531 -270.6674 -269.9366 -269.0702 -269.9689 -277.0594 -275.1290 -291.6407 -289.6269 -288.5471 -287.3730 -285.4722
 [3,] -275.16232 -277.2842 -275.9290 -275.3192 -274.3145 -270.6000 -269.6557 -267.7410 -270.6425 -270.8435 -268.8256 -289.6535 -289.5060 -287.1714 -289.0765 -283.6704
 [4,] -273.28643 -276.5743 -279.8084 -277.8297 -275.8461 -275.1641 -279.3690 -278.3020 -266.9166 -282.5261 -280.7962 -288.5563 -287.1536 -285.4661 -288.4532 -286.6260
 [5,] -273.72205 -274.7906 -280.4726 -275.8130 -281.6158 -275.5404 -277.9231 -281.6201 -269.8310 -281.9151 -280.3306 -286.8634 -284.5896 -283.1678 -285.0929 -282.7495
 [6,] -271.72639 -273.4516 -275.8477 -275.6839 -271.8633 -278.1729 -274.9926 -276.9595 -278.6623 -282.6819 -281.1888 -289.8140 -288.6669 -289.0610 -287.2544 -286.9919
 [7,] -271.70772 -271.8388 -278.5907 -279.6795 -280.1625        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA
 [8,]         NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA
 [9,]         NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA
[10,]         NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA
[11,]         NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA
[12,]         NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA
[13,]         NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA
[14,]         NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA
[15,]         NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA
[16,]         NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA


which.min(model.aic.test)  ## gives the index of the element in the matrix as if it were a long vector
# does not work here due to the NA values
## fortunatelly the following works:
which(model.aic==min(model.aic,na.rm=T),arr.ind=T)

##################################
###
## HOW do the ACF for AR and MA is supposed to look like?
## take r0=0; r_t ~ AR(4) coeff = 0.2 0.5 0.4 0.7. Is this stationary?

polynomial.rootsAR4=polyroot(c(1,0.2,0.5,0.1,0.07))
Mod(polynomial.rootsAR4)    ## we are fine

## use gaussian white noise with sd=0.3
## now to use the arima.sim function. the problem is that R uses the centered form.
## we need to recalculate the param: 

test.ar=arima.sim(n=365, model=list(order=c(4,0,0),ar=c(0.2,0.5,0.1,0.07),sd=0.3))
acf(test.ar)
pacf(test.ar)
ar(test.ar)
test.ma=arima.sim(n=365, model=list(order=c(0,0,4),ma=c(0.2,0.5,0.5,0.7),sd=0.3))
acf(test.ma)
pacf(test.ma)
acf(test.ma)$acf
test.arma=arima.sim(n=365, model=list(order=c(4,0,4),ar=c(0.2,0.5,0.1,0.07),ma=c(0.2,0.5,0.5,0.7),sd=0.3))
acf(test.arma)
pacf(test.arma)
ar(test.arma)

library(tseries) ## this is required by the function
eacf.tables(test.arma,max.p=8,max.q=8)  ## should give the corner at row 5, column 5 

## Another way: find the lowest AIC
max.p=8
max.q=8
model.aic=matrix(,nrow=max.p+1,ncol=max.q+1)
for (p in 0:max.p){
  for (q in 0:max.q){
  model.aic[p+1,q+1]= arima(test.arma,order=c(p,0,q))$aic
  }
}                   ## takes forever

print(model.aic)
which(model.aic==min(model.aic),arr.ind=T)