### Read in all city data
data_all <- read.csv(
  "http://www.stat.ufl.edu/~winner/data/weather_NOAA/allcities_annual_19602020.csv")
head(data_all)

### You will have 6 cities (These are read in alphabetical order) 
faa <- c("BOS","DFW","LAX","MIA", "ORD","SEA")

### Select the cities in your sample
city_sample <- data_all[data_all$city_FAA %in% faa,]  
head(city_sample); tail(city_sample)
attach(city_sample)

year_seq <- rep(1960:2020, times=6)  
city_seq <- rep(1:6, each=61)        

city_names <- c("Boston","Dallas","Los Angeles", 
                  "Miami","Chicago","Seattle")

## Make plots for cities 1-3 with ols (trend) lines added
# win.graph(height=5.5, width=7.0)
# quartz(height=5.5, width=7.0)

par(mfrow=c(2,3)) # Plot will have 2 rows and 3 columns 
# This strange sequence plots cities from East-West and North to South
for (i1 in c(6,5,1,3,2,4)) { 
  plot(aveTemp[city_seq == i1] ~ year_seq[city_seq == i1],
       main=paste("City = ", city_names[i1]),
       xlab="year", ylab="average annual temperature", col=(i1+1), pch=16)
  abline(lm(aveTemp[city_seq == i1] ~ year_seq[city_seq == i1]), 
    col=(i1+1), lwd=2)
}

par(mfrow=c(1,1))    # Reset to 1 plot per page

library(car)

mod.bos <- lm(aveTemp[city_seq==1] ~ I(year_seq[city_seq==1]-1960))
summary(mod.bos)

mod.dfw <- lm(aveTemp[city_seq==2] ~ I(year_seq[city_seq==2]-1960))
summary(mod.dfw)

mod.lax <- lm(aveTemp[city_seq==3] ~ I(year_seq[city_seq==3]-1960))
summary(mod.lax)

mod.mia <- lm(aveTemp[city_seq==4] ~ I(year_seq[city_seq==4]-1960))
summary(mod.mia)

mod.ord <- lm(aveTemp[city_seq==5] ~ I(year_seq[city_seq==5]-1960))
summary(mod.ord)

mod.sea <- lm(aveTemp[city_seq==6] ~ I(year_seq[city_seq==6]-1960))
summary(mod.sea)

durbinWatsonTest(mod.bos)
durbinWatsonTest(mod.dfw)
durbinWatsonTest(mod.lax)
durbinWatsonTest(mod.mia)
durbinWatsonTest(mod.ord)
durbinWatsonTest(mod.sea)

par(mfrow=c(2,3))
plot(resid(mod.sea), type="o", col=7, xlab="year", ylab="residual",
     main=paste("City = ", city_names[6]))
abline(h=0)
plot(resid(mod.ord), type="o", col=6, xlab="year", ylab="residual",
     main=paste("City = ", city_names[5]))
abline(h=0)
plot(resid(mod.bos), type="o", col=2, xlab="year", ylab="residual",
     main=paste("City = ", city_names[1]))
abline(h=0)
plot(resid(mod.lax), type="o", col=4, xlab="year", ylab="residual",
     main=paste("City = ", city_names[3]))
abline(h=0)
plot(resid(mod.dfw), type="o", col=3, xlab="year", ylab="residual",
     main=paste("City = ", city_names[2]))
abline(h=0)
plot(resid(mod.mia), type="o", col=5, xlab="year", ylab="residual",
     main=paste("City = ", city_names[4]))
abline(h=0)

par(mfrow=c(1,1))

detach(city_sample)
rm(list=ls(all=TRUE))