### Eksempler uke 9 require(astsa) ## Eks. 3.25 regr = ar.ols(rec, order=2, demean=FALSE, intercept=TRUE) fore = predict(regr, n.ahead=24) ts.plot(rec, fore$pred, col=1:2, xlim=c(1980,1990), ylab="Recruitment",ylim=c(0,120)) U = fore$pred+qnorm(0.975)*fore$se; L = pmax(fore$pred-qnorm(0.975)*fore$se,0) xx = c(time(U), rev(time(U))); yy = c(L, rev(U)) polygon(xx, yy, border = 8, col = gray(.6, alpha = .2)) lines(fore$pred, type="p", col=2) abline(h=mean(rec),lty=2,col=3) ## Eks. 3.26 x = arima.sim(list(order = c(1,0,1), ar =.9, ma=.5), n = 100) xr = rev(x) # xr is the reversed data pxr = predict(arima(xr, order=c(1,0,1)), 20) pxrp = rev(pxr$pred) pxrse = rev(pxr$se) nx = ts(c(pxrp, x), start=-19) plot(nx, ylab=expression(X[~t]), main='Backcasting') U = nx[1:20] + pxrse; L = nx[1:20] - pxrse xx = c(-19:0, 0:-19); yy = c(L, rev(U)) polygon(xx, yy, border = 8, col = gray(0.6, alpha = 0.2)) lines(-19:0, nx[1:20], col=2, type='o') abline(h=mean(x),lty=2,col=3) ## Eks. 3.28 rec.yw = ar.yw(rec, order=2) rec.yw$x.mean rec.yw$ar rec.yw$var.pred n = length(rec) muHat = mean(rec) rec0 = rec-muHat gammaHat0 = var(rec0)*(n-1)/n gammaHat1 = cov(rec0[-1],rec0[-n])*(n-2)/n gammaHat2 = cov(rec0[-(1:2)],rec0[-((n-1):n)])*(n-3)/n gammaVec = c(gammaHat1,gammaHat2) gammaMat = matrix(c(gammaHat0,gammaHat1,gammaHat1,gammaHat0),2,2) phiHat = solve(gammaMat,gammaVec) sigma2Hat = gammaHat0-t(gammaVec)%*%phiHat muHat phiHat sigma2Hat