mort = ts(scan("http://www.stat.ualberta.ca/~wiens/stat479/S&Sdatasets/cmort.dat"),start=1970, frequency=52) temp = ts(scan("http://www.stat.ualberta.ca/~wiens/stat479/S&Sdatasets/temp.dat"), start=1970, frequency=52) part = ts(scan("http://www.stat.ualberta.ca/~wiens/stat479/S&Sdatasets/part.dat"), start=1970, frequency=52) mort = ts(scan("cmort.dat"),start=1970, frequency=52) temp = ts(scan("temp.dat"), start=1970, frequency=52) part = ts(scan("part.dat"), start=1970, frequency=52) par(mfrow=c(3,1), mar=c(3,4,3,2)) plot(mort, main="Cardiovascular Mortality", xlab="", ylab="") plot(temp, main="Temperature", xlab="", ylab="") plot(part, main="Particulates", xlab="", ylab="") win.graph() # open another graphic device pairs(cbind(mort, temp, part)) # scatterplot matrix temp = temp-mean(temp) temp2 = temp^2 trend = 1:length(mort) # time win.graph() fit1 = lm(mort~ trend + temp + temp2 + part) summary(fit1) # regression results anova(lm(mort~cbind(trend, temp, temp2, part))) # pretty ANOVA table num = length(mort) # sample size AIC(fit1)/num - log(2*pi) # AIC (corresponds to def 2.1) AIC(fit1, k=log(num))/num - log(2*pi) # BIC (corresponds to def 2.3) (AICc = log(sum(resid(fit1)^2)/num)+(num+5)/(num-5-2)) # AICc (as in def 2.2) resids1 = fit1$residuals par(mfrow = c(2,1)) acf(resids1, xlab = "Lag (years)") pacf(resids1, xlab = "Lag (years)") residfit = arima(resids1, order= c(2,0,0), include.mean = F) phi1 = residfit$coef[1] phi1 # = .2184 phi2 = residfit$coef[2] phi2 # = .3623 # Transform: old_data_matrix = cbind(mort, trend, temp, temp2, part) new_data_matrix = filter(old_data_matrix, sides = 1, filter = c(1, -residfit$coef)) new_data_matrix = new_data_matrix[-c(1:length(residfit$coef)), ] new_mort = new_data_matrix[ ,1] new_trend = new_data_matrix[ ,2] new_temp = new_data_matrix[ ,3] new_temp2 = new_data_matrix[ ,4] new_part = new_data_matrix[ ,5] fit2 = lm(new_mort~ new_trend + new_temp + new_temp2 + new_part) summary(fit2) win.graph() resids2 = fit2$residuals par(mfrow = c(2,1)) acf(resids2, xlab = "Lag (weeks)") pacf(resids2, xlab = "Lag (weeks)") ############################ Try it once more? residfit2 = arima(resids2, order= c(1,0,0), include.mean = F) phi1 = residfit2$coef[1] # = 0.2184 new_data_matrix = filter(new_data_matrix, sides = 1, filter = c(1, -residfit2$coef)) new_data_matrix = new_data_matrix[-c(1:length(residfit2$coef)), ] new_mort = new_data_matrix[ ,1] new_trend = new_data_matrix[ ,2] new_temp = new_data_matrix[ ,3] new_temp2 = new_data_matrix[ ,4] new_part = new_data_matrix[ ,5] fit3 = lm(new_mort~ new_trend + new_temp + new_temp2 + new_part) summary(fit3) win.graph() resids3 = fit3$residuals par(mfrow = c(2,1)) acf(resids3, xlab = "Lag (weeks)") pacf(resids3, xlab = "Lag (weeks)")