1.3.1 Usual plot

plot(mytsdata, type="o", ylab="Quarterly Earnings")

1.3.2 ggplot

mytsdatafr<-data.frame(Y=as.matrix(mytsdata), date=as.Date(time(mytsdata)))
p <- ggplot(mytsdatafr, aes(x=date, y=Y)) +
  geom_line(alpha=1,  color="#69b3a2", size=2) +
  geom_point(alpha=1, color="red", size=2) +
  xlab("Time")+ylab("Quarterly Earnings")+
  scale_x_date(date_labels = "%Y %b %d")+
  theme(axis.text.x=element_text(angle=60, hjust=1))
ggplotly(p)

1.3.3 Plotly

data.frame(Y=as.matrix(mytsdata), date=as.Date(time(mytsdata))) %>%
plot_ly(x = ~date, y = ~Y, mode = 'lines+markers', type = "scatter", name="Confirm") %>% 
  layout(
    title = "",
    xaxis = list(rangeslider = list(type = "date")),
    yaxis = list(title = "Quarterly Earnings"))

1.3.4 dygraphs

library(dygraphs)
dygraph(mytsdata)%>% 
  dyOptions(drawPoints = TRUE, pointSize = 4)%>%
  dyHighlight(highlightSeriesOpts = list(strokeWidth = 3))%>% 
  dyRangeSelector()

1.4.1 Lag, Difference, and Autocorrelation

lag(mytsdata, 2)

           Jan       Feb       Mar       Apr       May       Jun       Jul
2009                                                                      
2010  1.000000 -1.839818 -3.017105 -1.842691  5.234162  9.841046  8.316776
2011 22.237127 16.637862  3.658659  7.686553 12.333570 17.485469 24.962183
2012 25.156549 26.243030 23.531001 20.708705 11.434616 12.531440 22.490719
2013 24.744799 24.425151 24.882243 36.510455 46.747396 47.942016 47.605089
2014 56.182945                                                            
           Aug       Sep       Oct       Nov       Dec
2009                                          0.000000
2010  6.511892  6.177134  8.671431 11.451401 19.952744
2011 28.333312 28.007451 30.626409 35.790274 38.104631
2012 37.753707 36.636157 36.766465 35.771098 25.376246
2013 56.433489 61.217564 48.260981 35.207093 40.722804
2014                                                  

diff(mytsdata, lag = 1, differences=1)

             Jan         Feb         Mar         Apr         May
2010                           1.0000000  -2.8398178  -1.1772875
2011   2.7799700   8.5013429   2.2843828  -5.5992647 -12.9792034
2012   5.1638647   2.3143578 -12.9480828   1.0864811  -2.7120285
2013  -0.9953678 -10.3948518  -0.6314469  -0.3196474   0.4570921
2014 -13.0538878   5.5157105  15.4601410                        
             Jun         Jul         Aug         Sep         Oct
2010   1.1744144   7.0768531   4.6068836  -1.5242693  -1.8048846
2011   4.0278940   4.6470171   5.1518996   7.4767133   3.3711296
2012  -2.8222960  -9.2740888   1.0968231   9.9592793  15.2629881
2013  11.6282120  10.2369407   1.1946202  -0.3369278   8.8284007
2014                                                            
             Nov         Dec
2010  -0.3347581   2.4942974
2011  -0.3258617   2.6189584
2012  -1.1175503   0.1303086
2013   4.7840752 -12.9565837
2014                        

acf(mytsdata, plot = FALSE)

Autocorrelations of series 'mytsdata', by lag

0.0000 0.0833 0.1667 0.2500 0.3333 0.4167 0.5000 0.5833 0.6667 0.7500 
 1.000  0.859  0.708  0.609  0.519  0.405  0.337  0.325  0.314  0.278 
0.8333 0.9167 1.0000 1.0833 1.1667 1.2500 1.3333 
 0.245  0.223  0.186  0.156  0.122  0.049 -0.016 

#acf(mytsdata, lag = 20)
#autoplot(acf(mytsdata, plot = FALSE), conf.int.fill = 'blue', conf.int.value = 0.5, conf.int.type = 'ma', colour="red")

1.4.2 Forcast Package

library(forecast)
ggtsdisplay(mytsdata)

gglagplot(mytsdata, do.lines = FALSE)

gglagchull(mytsdata)

ggAcf(mytsdata)

# Acf(mytsdata, lags = 9)
# acf1(mytsdata)
# lag.plot(mytsdata, lags = 9)

1.4.3 Trend

linearMod <- lm(mytsdata~time(mytsdata))
summary(linearMod)

Call:
lm(formula = mytsdata ~ time(mytsdata))

Residuals:
     Min       1Q   Median       3Q      Max 
-17.0685  -5.6420   0.6023   6.3084  16.8171 

Coefficients:
                  Estimate  Std. Error t value Pr(>|t|)    
(Intercept)    -23966.4357   1958.8331  -12.23 2.30e-16 ***
time(mytsdata)     11.9230      0.9735   12.25 2.22e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.278 on 48 degrees of freedom
Multiple R-squared:  0.7576,    Adjusted R-squared:  0.7525 
F-statistic:   150 on 1 and 48 DF,  p-value: 2.221e-16

anova(linearMod)

Analysis of Variance Table

Response: mytsdata
               Df  Sum Sq Mean Sq F value    Pr(>F)    
time(mytsdata)  1 10279.2 10279.2     150 2.221e-16 ***
Residuals      48  3289.4    68.5                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mytstrend <- linearMod$fitted.values
mytsde_trend <- linearMod$residuals

dataplot<-cbind(mytsdata, mytstrend, mytsde_trend) 
dygraph(dataplot)%>% 
  dyOptions(drawPoints = TRUE, pointSize = 4)%>%
  dyHighlight(highlightSeriesOpts = list(strokeWidth = 3))%>% 
  dyRangeSelector()

autoplot(linearMod, colour="red", which=1:6, ncol = 2, label.size = 3, title="Model checking plots")

shapiro.test(linearMod$residuals)

    Shapiro-Wilk normality test

data:  linearMod$residuals
W = 0.98465, p-value = 0.7566

1.4.4 Regression with Lagged Time Series

library(dynlm)
linearMod <- dynlm(mytsdata ~ L(mytsdata,1))
summary(linearMod)

Time series regression with "ts" data:
Start = 2010(3), End = 2014(3)

Call:
dynlm(formula = mytsdata ~ L(mytsdata, 1))

Residuals:
     Min       1Q   Median       3Q      Max 
-14.5796  -3.5823   0.0598   3.2617  15.4817 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)     2.72080    1.68965    1.61    0.114    
L(mytsdata, 1)  0.93266    0.05967   15.63   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.675 on 47 degrees of freedom
Multiple R-squared:  0.8387,    Adjusted R-squared:  0.8352 
F-statistic: 244.3 on 1 and 47 DF,  p-value: < 2.2e-16

dataplot<-cbind(mytsdata, trend=linearMod$fitted.values, residual=linearMod$residuals) 
dygraph(dataplot)%>% 
  dyOptions(drawPoints = TRUE, pointSize = 4)%>%
  dyHighlight(highlightSeriesOpts = list(strokeWidth = 3))%>% 
  dyRangeSelector()

1.4.5 Seasonality

ggtsdisplay(mytsde_trend)

ggtsdisplay(diff(mytsde_trend, lag = 1, differences=6))

library(forecast)
ts.stl <- stl(mytsdata, "periodic")  # decompose the TS
ts.sa  <- seasadj(ts.stl)
autoplot(ts.stl) 

dataplot<-cbind(mytsdata, seasonal=ts.stl$time.series[,1], trend=ts.stl$time.series[,2]) 
dygraph(dataplot)%>% 
  dyOptions(drawPoints = TRUE, pointSize = 4)%>%
  dyHighlight(highlightSeriesOpts = list(strokeWidth = 3))%>% 
  dyRangeSelector()

seasonplot(ts.sa, 5, col=rainbow(5), year.labels=TRUE, main="Seasonal plot: Airpassengers", lwd=3)

1.4.6 Making Stationary

library(tseries)
adf.test(mytsdata) # p-value < 0.05 indicates the TS is stationary

    Augmented Dickey-Fuller Test

data:  mytsdata
Dickey-Fuller = -3.3154, Lag order = 3, p-value = 0.07932
alternative hypothesis: stationary

nsdiffs(mytsdata)  # number for seasonal differencing needed

[1] 0

ndiffs(mytsdata) # number of differences need to make it stationary

[1] 1

stationaryTS <- diff(mytsdata, differences= ndiffs(mytsdata))
adf.test(stationaryTS) # p-value < 0.05 indicates the TS is stationary

    Augmented Dickey-Fuller Test

data:  stationaryTS
Dickey-Fuller = -3.9373, Lag order = 3, p-value = 0.01973
alternative hypothesis: stationary

dataplot<-cbind(mytsdata,stationaryTS) 
dygraph(dataplot)%>% 
  dyOptions(drawPoints = TRUE, pointSize = 4)%>%
  dyHighlight(highlightSeriesOpts = list(strokeWidth = 3))%>% 
  dyRangeSelector()

1.4.7 Forecast

ts.stl %>% forecast() %>%
  autoplot()

fit <- auto.arima(ts.stl$time.series[,3], seasonal=FALSE)
summary(fit)

Series: ts.stl$time.series[, 3] 
ARIMA(2,0,0) with zero mean 

Coefficients:
         ar1      ar2
      0.8735  -0.6074
s.e.  0.1167   0.1158

sigma^2 estimated as 18.39:  log likelihood=-143.36
AIC=292.72   AICc=293.24   BIC=298.45

Training set error measures:
                    ME     RMSE      MAE      MPE    MAPE      MASE
Training set 0.1290637 4.201792 3.438581 124.7528 161.216 0.4287052
                   ACF1
Training set 0.03800343

fit %>% forecast() %>%
  autoplot()

1.4.8 Entir Process

mytsdata %>% stl(s.window='periodic') %>% seasadj() -> mytsdatadj
autoplot(mytsdatadj)

The data are clearly non-stationary, as the series wanders up and down for long periods. Consequently, we will take a first difference of the data.

mytsdatadj %>% diff() %>% ggtsdisplay(main="")

The PACF shown is suggestive of an AR(2) model. So an initial candidate model is an ARIMA(2,1,0). There are no other obvious candidate models.

(fit <- Arima(mytsdatadj, order=c(2,1,0)))

Series: mytsdatadj 
ARIMA(2,1,0) 

Coefficients:
         ar1      ar2
      0.3716  -0.4716
s.e.  0.1341   0.1319

sigma^2 estimated as 28.68:  log likelihood=-151.02
AIC=308.04   AICc=308.57   BIC=313.71

checkresiduals(fit)

    Ljung-Box test

data:  Residuals from ARIMA(2,1,0)
Q* = 7.5718, df = 8, p-value = 0.4764

Model df: 2.   Total lags used: 10

autoplot(forecast(fit))

1.4.9 Seasonal Package

library(seasonal)
m <- seas(mytsdata)
#view(m)

dataplot<-cbind(trend=m[["data"]][,"trend"], seasonaladj=m[["data"]][,"seasonaladj"], irregular=m[["data"]][,"irregular"],  data=mytsdata)
  dygraph(dataplot)%>%
    dyOptions(gridLineColor = "#E1E5EA", axisLineColor = "#303030", drawPoints
=TRUE, strokeWidth = 2, colors = c("blue","green","grey","red"))%>%
    dyHighlight(highlightSeriesOpts = list(strokeWidth = 3))%>% 
    dyRangeSelector()

1.5.1 White Noise

n <- 200
value <-  rnorm(n,0,1)
mytsdata <- ts(value, start = 1, end = n)
mytsdatafr<-data.frame(Y=as.matrix(mytsdata), date=as.Date(time(mytsdata)))
p <- ggplot(mytsdatafr, aes(x=date, y=Y)) +
  geom_line(alpha=1,  color="#69b3a2", size=2) +
  geom_point(alpha=1, color="red", size=2) +
  xlab("Time")+ylab("Value")
ggplotly(p)

autoplot(acf(mytsdata, plot = FALSE))

Smoothing

values <- stats::filter(value, sides=2, filter=rep(1/5,5))
mytsdata <- ts(values, start = 1, end = n)
mytsdatafr<-data.frame(Y=as.matrix(mytsdata), date=as.Date(time(mytsdata)))
p <- ggplot(mytsdatafr, aes(x=date, y=Y)) +
  geom_line(alpha=1,  color="#69b3a2", size=2) +
  geom_point(alpha=1, color="red", size=2) +
  xlab("Time")+ylab("Value")
ggplotly(p)

autoplot(acf(mytsdata[-c(1,2,199,200)], plot = FALSE))

1.5.2 Trend and White Noise

value <- (1:n)/10 + rnorm(n,0,1)
mytsdata <- ts(value, start = 1, end = n)
mytsdatafr<-data.frame(Y=as.matrix(mytsdata), date=as.Date(time(mytsdata)))
p <- ggplot(mytsdatafr, aes(x=date, y=Y)) +
  geom_line(alpha=1,  color="#69b3a2", size=2) +
  geom_point(alpha=1, color="red", size=2) +
  xlab("Time")+ylab("Value")
ggplotly(p)

autoplot(acf(mytsdata, plot = FALSE))

1.5.3 Autoregresive

value <- 0
for (i in 2:n) {
  value[i] <- value[i-1] + rnorm(1,0,1)
}
mytsdata <- ts(value, start = 1, end = n)
mytsdatafr<-data.frame(Y=as.matrix(mytsdata), date=as.Date(time(mytsdata)))
p <- ggplot(mytsdatafr, aes(x=date, y=Y)) +
  geom_line(alpha=1,  color="#69b3a2", size=2) +
  geom_point(alpha=1, color="red", size=2) +
  xlab("Time")+ylab("Value")
ggplotly(p)

autoplot(acf(mytsdata, plot = FALSE))

value <- 0
value[2] <- 1
for (i in 3:n) {
  value[i] <- value[i-1] - value[i-2]/9 + rnorm(1,0,1)
}
mytsdata <- ts(value, start = 1, end = n)
mytsdatafr<-data.frame(Y=as.matrix(mytsdata), date=as.Date(time(mytsdata)))
p <- ggplot(mytsdatafr, aes(x=date, y=Y)) +
  geom_line(alpha=1,  color="#69b3a2", size=2) +
  geom_point(alpha=1, color="red", size=2) +
  xlab("Time")+ylab("Value")
ggplotly(p)

autoplot(acf(mytsdata, plot = FALSE))

1.5.4 Seasonality and White Noise

value <- 2*cos(2*pi*1:n/10 +0.6*pi) + rnorm(n,0,1)
mytsdata <- ts(value, start = 1, end = n)
mytsdatafr<-data.frame(Y=as.matrix(mytsdata), date=as.Date(time(mytsdata)))
p <- ggplot(mytsdatafr, aes(x=date, y=Y)) +
  geom_line(alpha=1,  color="#69b3a2", size=2) +
  geom_point(alpha=1, color="red", size=2) +
  xlab("Time")+ylab("Value")
ggplotly(p)

autoplot(acf(mytsdata, plot = FALSE))

1.5.5 Cross-correlation

value <- (1:n)/20 + cos((1:n)/10) + rnorm(n,0,1)
mytsdatax <- ts(value, start = 1, end = n)
value <- (1:n)/20 + rnorm(n,0,2) + 5
mytsdatay <- ts(value, start = 1, end = n)
mytsdata <- cbind(mytsdatax, mytsdatay)
dygraph(mytsdata)%>% 
  dyOptions(drawPoints = TRUE, pointSize = 4)%>%
  dyHighlight(highlightSeriesOpts = list(strokeWidth = 3))%>% 
  dyRangeSelector()

ggCcf(mytsdatax, mytsdatay)

library(astsa) # Applied Statistical Time Series Analysis

Attaching package: 'astsa'

The following object is masked from 'package:seasonal':

    unemp

The following object is masked from 'package:forecast':

    gas

lag2.plot (mytsdatax, mytsdatay, 11)