data("UKgas")
? UKgas

ts.plot(UKgas,type="o",main="UKgas")
ts.plot(log(UKgas),type="o",main="log UKgas")

# log trend model
y = log(UKgas)
n = length(UKgas)
trend = 1:n
m1 = lm(y~trend)
summary(m1)$coef
summary(m1)$r.squared

# add dummies for quarters
quarter = rep(1:4,27)
D2 = as.numeric(quarter==2)
D3 = as.numeric(quarter==3)
D4 = as.numeric(quarter==4)
summary(lm(y~trend+D2+D3+D4))$coef

# A shortcut
m2 = lm(y~trend+factor(quarter))
summary(m2)$coef
summary(m2)$r.squared

predict(m2,data.frame(trend=(n+1),quarter=1))
predict(m2,data.frame(trend=(n+2),quarter=2))
co = coef(m2)
co[1]+co[2]*(n+1)
co[1]+co[2]*(n+2) + co[3]

# add lag value of y
ylag1 = c(NA, y[1:(n-1)])
m3 = lm(y~trend+factor(quarter)+ylag1)
summary(m3)$coef
summary(m3)$r.squared

yhat = predict(m3,data.frame(trend=(n+1),quarter=1,ylag1=y[n]))
yhat
yhat = predict(m3,data.frame(trend=(n+2),quarter=2,ylag1=yhat))
yhat

# try fourth lag
ylag4 = c(NA, NA, NA, NA, y[1:(n-4)])
m4 = lm(y~trend+factor(quarter)+ylag4)
summary(m4)$coef
summary(m4)$r.squared

predict(m4,data.frame(trend=(n+1),quarter=1,ylag4=y[n-3]))
predict(m4,data.frame(trend=(n+2),quarter=2,ylag4=y[n-2]))

# optional advanced stuff
# predict ukgas based on model for log ukgas
exp(7.128523+summary(m4)$sigma^2/2)
exp(6.515706+summary(m4)$sigma^2/2)

# ar(1) for seasonal difference
dy4 = diff(y,4)
dy4lag = c(NA,dy4[1:(length(dy4)-1)])
summary(lm(dy4~dy4lag))$coef
plot(dy4,main="seasonal difference of log UKgas")

# acf---detecting seasonality and helping model specification
acf(UKgas)
acf(diff(UKgas))
acf(diff(UKgas,4))