## 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.08=read.csv("GSPCmonthly2000-June2008.csv",header=T) sp500monthly00.08.price= sp500monthly00.08$Adj.Close[length(sp500monthly00.08$Adj.Close):1] ## 102 months library(Rmetrics) #TO calculate returns automatically: returns( sp500monthly00.08.price) ## default gives continuously compunded returns. Note the first obs is NA returns( sp500monthly00.08.price,method="compound") ## same as above returns( sp500monthly00.08.price,method="simple") ## gives the simple returns returns( sp500monthly00.08.price,method="discrete") ## same as simple returns above sp500monthly00.08.return= returns( sp500monthly00.08.price)[-1] ## to eliminate the first NA element ## To calculate overlapping returns look at documentation: ?RollingAnalysis # we need to apply sum function: sp500overlapping.returns= rollFun(sp500monthly00.08.return,n=12,FUN=sum) acf(sp500overlapping.returns) ## Note the MA(11) 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 print(model.aic) Error in optim(init[mask], armafn, method = "BFGS", hessian = TRUE, control = optim.control, : non-finite finite-difference value [11] In addition: There were 40 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,] -76.70684 -165.9827 -195.5093 -227.2617 -236.5673 -241.8049 -242.9804 -248.7442 -252.0641 -257.9901 -257.0120 -280.5612 -279.4618 -277.4635 -276.5734 -275.2790 [2,] -259.42225 -257.9995 -256.0232 -255.8078 -254.8367 -253.5283 -255.4713 -252.7425 -255.6870 -264.6866 -262.8027 -279.4611 -277.7457 -275.8822 -274.5686 -272.8864 [3,] -257.99477 -260.7203 -258.7234 -253.6190 -255.4398 -262.1020 -260.8965 -248.8004 -260.9134 -263.3517 -261.2039 -277.4613 -275.4614 -278.8392 -273.5066 -271.5256 [4,] -255.99679 -258.7240 -258.3112 -261.8293 -258.1397 -256.5045 -258.3013 -260.9694 -259.1329 -262.4585 -261.0639 -275.7691 -274.1143 -273.7186 -276.1807 -274.8691 [5,] -255.51270 -254.9874 -261.9078 -260.3628 -258.8598 -256.4526 -255.3503 -264.2987 -258.9272 -263.7503 -256.7379 -274.2170 -275.1103 -276.6504 -271.0130 -273.1628 [6,] -254.72329 -255.5516 -261.3770 -260.7642 -257.6201 -257.6089 -259.7535 -258.2091 -262.1222 -268.7683 -255.9120 -279.2880 -277.2960 -275.7860 -274.0940 -273.6778 [7,] -254.54885 -253.7864 -257.5772 -249.3308 -255.7179 -253.7004 -258.8523 -260.9512 -257.6784 -267.7631 -266.0649 -277.3001 -275.6032 -273.8687 -271.6268 -269.7126 [8,] -253.28742 -251.6702 -251.9445 -255.6907 -257.9198 -254.3611 -256.8495 -257.9765 -260.3036 -259.8824 -258.6235 -276.0570 -273.8345 -272.4767 -271.7395 -272.6487 [9,] -251.33219 -249.3752 -250.3624 -253.7650 -255.3460 -251.1157 -256.2748 -264.0831 -261.7783 -257.9613 -256.3976 -274.2211 -272.0788 -270.5796 -270.8024 -270.8564 [10,] -250.85347 -248.9227 -251.5438 -251.9959 -261.3308 -248.1957 -253.4628 -259.4627 -260.1809 -268.7085 -267.2504 -272.7133 -270.6654 -268.8124 -269.6696 -269.3095 [11,] -248.91969 -251.8737 -254.7337 -256.3304 -259.7102 -251.2660 -252.9392 -262.3321 -261.3983 -266.9446 -260.2334 -270.8820 -268.7318 -266.8471 -268.0403 -265.9874 [12,] -247.03602 -245.0567 -252.7403 -254.1501 -256.2900 -261.5867 -262.2474 -261.6258 -257.5392 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) ## gives the index of the element in the matrix as if it were a long vector which(model.aic==min(model.aic),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)