Apply the Box-Jenkins method (transform to stationarity if necessary and identify the time series model(e.g., ARMA(p,q))) for the luteinizing hormone data. Write down the algebraic expression of the fitted model. How is your "final" model compared with an AR(3) ? Use your final model to forecast the next 24 luteinizing hormone measurements

Code EXAMPLE FROM BOOK

######################################################## #3.5: ESTIMATION/FITTING ########################################################

library(astsa) #Example 3.27(text). AR(2): xt=1.5x(t-1)-.75x(t-2) +wt set.seed(8675309) ar2.sim=arima.sim(list(order = c(2,0,0), ar = c(1.5,-0.75)), n = 144) acf2(ar2.sim)

########### #estimation with goodness-of-fit diagnostics sarima(ar2.sim, 2,0,0) sarima(ar2.sim, 2,0,0, no.constant=TRUE) #Example 3.29 (text): MA(1), xt= wt + theta w(t-1) , theta=0.9. set.seed(2) ma1 = arima.sim(list(order = c(0,0,1), ma = 0.9), n = 150)

########### #estimation with goodness-of-fit diagnostics sarima(lh, 1,0,0) sarima(lh, 1,0,0, no.constant=TRUE)

#MODEL: x[t] = 0.8413*w[t-1] +w[t]

#ARMA(1,1) , phi =0.83, theta= -0.43. set.seed(2786) arma11.sim = arima.sim(list(order = c(1,0,1), ar=0.83, ma = -0.43), n = 150) ########### #estimation with goodness-of-fit diagnostics sarima(arma11.sim,1,0,0) #MODEL: Subtract muhat from x[t] only like #(x[t]- (-0.2487) ) = 0.7958*(x[t-1]-(-0.2487)) - 0.2913*w[t-1] + w[t]

#But since the p-value is large for 'xmean' (meaning mu_xt = 0) then sarima(arma11.sim, 1,0,0, no.constant=TRUE)

#Hence, MODEL: x[t] = 0.7958*x[t-1] - 0.2913*w[t-1] + w[t]

############# #Recruitment data #estimation with goodness-of-fit diagnostics acf2(rec) sarima(rec, 13,0,0) #Hence, MODEL: x[t]-61.85 = 1.35*(x[t-1]-61.85) - 0.46*(x[t-2]-61.85) + w[t] #OR sarima(rec, 13,0,0, fixed=c(NA,NA,0,0,0,0,0,0,0,0,0,0,NA, NA) )

#OR sarima(rec, 13,0,0, fixed=c(NA,NA,0,NA,NA,0,0,0,0,0,0,0,NA, NA) )

#auto.arima auto.arima(rec) sarima(rec,1,1,0,0,0,2,12, no.constant=TRUE)

#Monthly lung disease deaths for females (fdeaths) #auto.arima acf2(fdeaths) sarima(fdeaths, 4,0,0 ) #OR sarima(rec, 4,0,0, fixed=c(NA,0,0,NA, NA) )

auto.arima(fdeaths) sarima(fdeaths,0,0,0,2,1,0,12, no.constant=FALSE)

########################################### #Section 3.4: FORECASTING # ###########################################

#Global warming data: temp deviations (1880-2015) library(astsa) data(globtemp) plot(globtemp) str(globtemp) time=time(globtemp) time2=time^2 time(globtemp) dat=cbind(globtemp,time,time2) mod1=lm(globtemp~I(time)+ I(time2), dat) mod1 #plot fit plot(time,globtemp, xlim=c(1880,2015), ylim=c(-0.5, 1),frame.plot=FALSE) lines( 2.844e+02 -2.992e-01*time +7.860e-05*time2, lwd=3, col="red")

#residuals qqnorm(resid(mod1),frame.plot=FALSE) qqline(resid(mod1))

#prediction for 2050 predict(mod1, data.frame(time=2050, time2=2050^2 ) , interval = "conf",level=0.95)

#Recruitment data library(astsa) sarima.for(rec,24,2,0,0) #with +/- 1,2*SE (prediction error bounds)

#OR library(astsa) sarima.for(fdeaths,24,0,0,0,2,1,0,12) #with +/- 1,2*SE (prediction error bounds)