TS in R

##Decomposition ## TS Class

AP <- AirPassengers
AP
##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1949 112 118 132 129 121 135 148 148 136 119 104 118
## 1950 115 126 141 135 125 149 170 170 158 133 114 140
## 1951 145 150 178 163 172 178 199 199 184 162 146 166
## 1952 171 180 193 181 183 218 230 242 209 191 172 194
## 1953 196 196 236 235 229 243 264 272 237 211 180 201
## 1954 204 188 235 227 234 264 302 293 259 229 203 229
## 1955 242 233 267 269 270 315 364 347 312 274 237 278
## 1956 284 277 317 313 318 374 413 405 355 306 271 306
## 1957 315 301 356 348 355 422 465 467 404 347 305 336
## 1958 340 318 362 348 363 435 491 505 404 359 310 337
## 1959 360 342 406 396 420 472 548 559 463 407 362 405
## 1960 417 391 419 461 472 535 622 606 508 461 390 432
?AirPassengers
## starting httpd help server ... done
cl <- class(AP)
st <- start(AP)
e <- end(AP)
fr <- frequency(AP)
cl
## [1] "ts"
st
## [1] 1949    1
e
## [1] 1960   12
fr
## [1] 12
plot(AP, ylab = "Passengers (1000's)")

layout(matrix(c(1,2), ncol=1))
plot(aggregate(AP), ylab = "Aggregated annually")
boxplot(AP ~ cycle(AP))

Maine.month <- read.table("Maine.dat", header=TRUE)
attach(Maine.month)
class(Maine.month)
## [1] "data.frame"
head(Maine.month)
##   unemploy
## 1      6.7
## 2      6.7
## 3      6.4
## 4      5.9
## 5      5.2
## 6      4.8
layout(matrix(c(1,2), ncol=1))
Maine.month.ts <- ts(unemploy, start = c(1996,1), freq = 12)
plot(Maine.month.ts, ylab = "unemployed (%)", main = " Monthly:January 1996-August 2006")
Maine.annual <- aggregate(Maine.month.ts)
plot(Maine.annual, ylab = "unemployed(%)", main = "Annual: 1996-2006")

Maine.Feb <- window(Maine.month.ts, start = c(1996,2), freq = TRUE)
Maine.Aug <- window(Maine.month.ts, start = c(1996,8), freq = TRUE)
Feb.ratio <- mean(Maine.Feb) / mean(Maine.month.ts)
Aug.ratio <- mean(Maine.Aug) / mean(Maine.month.ts)
Feb.ratio
## [1] 1.222529
Aug.ratio
## [1] 0.8163732

###Multiple time series

cbe <- read.table("cbe.dat", header=TRUE)
attach(cbe)
class(cbe)
## [1] "data.frame"
head(cbe)
##   choc beer elec
## 1 1451 96.3 1497
## 2 2037 84.4 1463
## 3 2477 91.2 1648
## 4 2785 81.9 1595
## 5 2994 80.5 1777
## 6 2681 70.4 1824
Choc.ts <- ts(cbe$choc, start = c(1958,1), frequency = 12)
Beer.ts <- ts(cbe$beer, start = c(1958,1), frequency = 12)
Elec.ts <- ts(cbe$elec, start = c(1958,1), frequency = 12)

combined.ts <- cbind(Elec.ts, Beer.ts, Choc.ts)

plot(combined.ts, main = "cbind(Elec.ts, Beer.ts, Choc.ts)")

AP.elec <- ts.intersect(AP, Elec.ts)
AP <- AP.elec[,1]
Elec <- AP.elec[,2]
layout(1:3) 
  
plot(AP, main = "", ylab = "Air passengers / 1000's") 
plot(Elec, main = "", ylab = "Electricity production / MkWh")
plot(as.vector(AP), as.vector(Elec), xlab = "Air passengers / 1000's", ylab = "Electricity production / MWh")
abline(reg = lm(Elec ~ AP))

# plot(decompose(Elec.ts))
Elec.decom <- decompose(Elec.ts, type = "additive")
plot(Elec.decom)

Trend <- Elec.decom$trend
Seasonal <- Elec.decom$seasonal
ts.plot(cbind(Trend, Trend + Seasonal), lty = 1:2)

plot(stl(Elec.ts, s.window = 7, t.jump = 1))

plot(stl(Elec.ts, s.window =7))

#plot(stl(Elec.ts,"per"))
plot(stl(Elec.ts, s.window = "periodic"))

##Regression

library(forecast)
## Warning: package 'forecast' was built under R version 4.5.1
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(ggplot2)

?mdeaths
class(mdeaths)
## [1] "ts"
start(mdeaths)
## [1] 1974    1
end(mdeaths)
## [1] 1979   12
frequency(mdeaths)
## [1] 12
# autoplot  of a ts object
autoplot(mdeaths)

# autoplot of a forecase object
fc <- forecast(fdeaths)
autoplot(fc)

# autoplot of an stl object
autoplot(stl(mdeaths, s.window = "periodic", robust = TRUE))

# plotting multiple forecasts in on plot
fmdeaths <- cbind(Males = mdeaths, Females = fdeaths)
fit <- tslm(fmdeaths ~ trend + season)
fcast <- forecast(fit, h = 10)
autoplot(fcast)

# plotting the components of an ETS model
fit <- ets(mdeaths)
autoplot(fit)

# plotting the inverse characteristic roots of an ARIMA model
fit <- auto.arima(mdeaths, D=1)
autoplot(fit)

ggtsdisplay(mdeaths)

ggAcf(mdeaths)

ggPacf(mdeaths)

ggCcf(mdeaths, fdeaths)

ggseasonplot(mdeaths)

ggmonthplot(mdeaths)

autoplot(mdeaths) + geom_forecast(h=36)

# TSLM
deaths.lm <- tslm(mdeaths ~ trend + season)
summary(deaths.lm)
## 
## Call:
## tslm(formula = mdeaths ~ trend + season)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -422.70  -58.60   -0.11   64.29  644.03 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2257.1179    74.8463  30.157  < 2e-16 ***
## trend          -4.1060     0.9662  -4.250 7.72e-05 ***
## season2       -44.3940    97.0086  -0.458  0.64890    
## season3      -151.1214    97.0230  -1.558  0.12468    
## season4      -460.1821    97.0471  -4.742 1.38e-05 ***
## season5      -799.2429    97.0807  -8.233 2.21e-11 ***
## season6      -922.4702    97.1240  -9.498 1.71e-13 ***
## season7      -968.5310    97.1768  -9.967 2.91e-14 ***
## season8     -1063.4250    97.2393 -10.936 8.09e-16 ***
## season9     -1075.3190    97.3112 -11.050 5.34e-16 ***
## season10     -851.7131    97.3927  -8.745 3.04e-12 ***
## season11     -711.1071    97.4838  -7.295 8.57e-10 ***
## season12     -288.1679    97.5843  -2.953  0.00451 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 168 on 59 degrees of freedom
## Multiple R-squared:  0.875,  Adjusted R-squared:  0.8495 
## F-statistic: 34.41 on 12 and 59 DF,  p-value: < 2.2e-16
# resampling
quarterly <- aggregate(mdeaths, nfrequency = 4)
deathsq.lm <- tslm(quarterly ~ trend + season)
summary(deathsq.lm)
## 
## Call:
## tslm(formula = quarterly ~ trend + season)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -677.20 -145.69   43.33  144.09  797.43 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6588.16     175.57  37.524  < 2e-16 ***
## trend         -36.95      10.01  -3.692  0.00155 ** 
## season2     -1986.38     193.66 -10.257 3.49e-09 ***
## season3     -2911.76     194.43 -14.976 5.66e-12 ***
## season4     -1655.47     195.72  -8.458 7.25e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 335 on 19 degrees of freedom
## Multiple R-squared:  0.9325, Adjusted R-squared:  0.9183 
## F-statistic: 65.59 on 4 and 19 DF,  p-value: 7.482e-11
autoplot(stl(quarterly, "per"))

# Forecasting
data = AirPassengers
# creat samples
training = window(data, start = c(1949,1), end = c(1955,12))
validation = window(data, start = c(1956,1))
# naive method
library(MLmetrics)
## Warning: package 'MLmetrics' was built under R version 4.5.1
## 
## Attaching package: 'MLmetrics'
## The following object is masked from 'package:base':
## 
##     Recall
naive = naive(training, h = length(validation))
nmap = MAPE(naive$mean, validation) * 100
nmap
## [1] 27.27412
snaive = snaive(training, h = length(validation))
smap = MAPE(snaive$mean, validation) * 100
smap
## [1] 27.04689
snaive_fit <- snaive(training, h = length(validation))
plot(data, col = "blue", xlab = "Year", ylab = "Passengers", main = "Seasonal Naive Forecast", type='l')
lines(snaive_fit$mean, col = "red", lwd=2)

# arima model
arima_optimal = auto.arima(training)
print(arima_optimal)
## Series: training 
## ARIMA(0,1,1)(1,1,0)[12] 
## 
## Coefficients:
##           ma1     sar1
##       -0.2591  -0.2603
## s.e.   0.1296   0.1179
## 
## sigma^2 = 96.77:  log likelihood = -262.5
## AIC=531   AICc=531.36   BIC=537.79
library(astsa)
## Warning: package 'astsa' was built under R version 4.5.1
## 
## Attaching package: 'astsa'
## The following object is masked from 'package:forecast':
## 
##     gas
sarima_forecast = sarima.for(training, n.ahead=length(validation),p=0,d=1,q=1,P=1,D=1,Q=0,S=12)

MAPE(sarima_forecast$pred, validation) * 100
## [1] 6.544624