Chapter 3 Time Series Data Pre-Processing and Visualization

rm(list = ls())
setwd("C:/Users/Tejendra/Desktop/FoldersOnDesktop/UdemyCourse/Section3")
require(tidyverse)
require(tidymodels)
require(data.table)
require(tidyposterior)
require(tsibble)  #tsibble for time series based on tidy principles
require(fable)  #for forecasting based on tidy principles
require(ggfortify)  #for plotting timeseries
require(forecast)  #for forecast function
require(tseries)
require(chron)
require(lubridate)
require(directlabels)
require(zoo)

setwd("C:/Users/Tejendra/Desktop/FoldersOnDesktop/UdemyCourse/Section3")
## U Standard time series functions
# Lets create a time series object - class ts
# Getting data
mydata = runif(n = 50, min = 10, max = 45)

# ts for class time series
# Data starts in 1956 - 4 observations/year (quarterly)
mytimeseries = ts(data = mydata, 
                  start = 1956, frequency = 4)

# Lets see how the data looks
plot(mytimeseries)

# Checking the class
class(mytimeseries)

## [1] "ts"

# Checking the timestamp
time(mytimeseries)

##         Qtr1    Qtr2    Qtr3    Qtr4
## 1956 1956.00 1956.25 1956.50 1956.75
## 1957 1957.00 1957.25 1957.50 1957.75
## 1958 1958.00 1958.25 1958.50 1958.75
## 1959 1959.00 1959.25 1959.50 1959.75
## 1960 1960.00 1960.25 1960.50 1960.75
## 1961 1961.00 1961.25 1961.50 1961.75
## 1962 1962.00 1962.25 1962.50 1962.75
## 1963 1963.00 1963.25 1963.50 1963.75
## 1964 1964.00 1964.25 1964.50 1964.75
## 1965 1965.00 1965.25 1965.50 1965.75
## 1966 1966.00 1966.25 1966.50 1966.75
## 1967 1967.00 1967.25 1967.50 1967.75
## 1968 1968.00 1968.25

# Refining the start argument
mytimeseries = ts(data = mydata, 
                  start = c(1956,3), frequency = 4)
## Creating a ts object - Exercise
# Get a random walk of 450 numbers, eg rnorm, runif, etc 
# In the solution I am going to use a cumulative sum on the normal distribution x = cumsum(rnorm(n = 450)) 
# If you want it to be reproducible, you can set a seed 
# Add the time component: it is a monthly dataset, which starts in November 1914 
# Get a simple plot for this time series # Advanced: how would you get the same type with the "lattice" package?

set.seed(2019)
ts1 <- cumsum(rnorm(n = 450))
ts1 <- ts(ts1,
          start = c(1914, 11),
          frequency = 12)
plot(ts1)

lattice::xyplot.ts(ts1)

## U Plots for time series data
# Standard R Base plots
plot(nottem) 

# Plot of components
plot(decompose(nottem)) 

# Directly plotting a forecast of a model
plot(forecast(auto.arima(nottem)), h = 5)

## Warning in plot.window(xlim, ylim, log, ...): "h" is not a graphical
## parameter

## Warning in title(main = main, xlab = xlab, ylab = ylab, ...): "h" is not a
## graphical parameter

## Warning in box(...): "h" is not a graphical parameter

# Random walk
plot.ts(cumsum(rnorm(500)))

# The ggplot equivalent to plot
autoplot((nottem))

# Ggplots work with different layers
autoplot(nottem) + ggtitle("Autoplot of Nottingham temperature data") +
  theme_minimal()

# Time series specific plots
ggseasonplot(nottem) + theme_minimal()

ggmonthplot(nottem)

## Exercise Seasonplot - library (forecast)
# use the seasonplot function in order  
# to resemble the plot shown here 
# we are using the AirPassengers dataset  
# for all the needed arguments on the plot, 
# check out the help for par 
# make sure that the labels are visible

AirPassengers

##      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

seasonplot(AirPassengers)

seasonplot(AirPassengers, xlab = "", col = c("red", "blue"), year.labels = T, labelgap = 0.35, type = "s",  bty = "7", cex = 0.75,  main = "Seasonal plot of dataset AirPassengers")

theme_set(theme_minimal())
ggseasonplot(AirPassengers, 
             year.labels = TRUE) + 
  xlab("") +
  ggtitle("Seasonal plot of dataset AirPassengers") +
  theme(plot.title = element_text(hjust = 0.5)) + 
  scale_color_manual(values=rep(c("red", "blue"), 6)) 

## German Inflation Rates
#https://www.statbureau.org/en/germany/inflation-tables

# data is copied from web, first and last col deleted
#mydata = scan()
#plot.ts(mydata)
#germaninfl = ts(mydata, start = 2008, frequency = 12)
#plot(germaninfl)

### Working with Irregular Time Series
## dataset: irregular_sensor

## Method 1 - removing the time component
irregular_sensor <- read_csv("./Data/irregular-sensor.csv", 
                             col_names = FALSE)

## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   X2 = col_double()
## )

# Irregular_sensor time series csv
class(irregular_sensor$X1) 

## [1] "character"

#separate a column into two or more components
#using separate function from tidyr
irreg.split = tidyr::separate(irregular_sensor, 
                              col = X1, #column to separate
                       into = c('date', 'time'), #new column names
                       sep = 8, # where to separate from
                       remove = T) #remove the original column after separation

# Using only the date
sensor.date = strptime(irreg.split$date, 
                       '%m/%d/%y') #format of the sepapration 

# Creating a data.frame for orientation
irregts.df = data.frame(date = as.Date(sensor.date), 
                        measurement = irregular_sensor$X2)
irregts.df = data_frame(1:25) %>% 
  mutate(date = as.Date(sensor.date),
         measurement = irregular_sensor$X2) %>% 
  select(-`1:25`) %>% 
  group_by(date) %>% 
  mutate(measurement = mean(measurement)) %>% 
  arrange(date) %>% 
  distinct()

## Error in select(., -`1:25`): unused argument (-`1:25`)

# Getting a zoo object
irreg.dates = zoo(irregts.df$measurement,
                  order.by = irregts.df$date)

## Warning in zoo(irregts.df$measurement, order.by = irregts.df$date): some
## methods for "zoo" objects do not work if the index entries in 'order.by'
## are not unique

# Regularizing with aggregate
ag.irregtime = aggregate(irreg.dates, 
                         as.Date, mean)
ag.irregtime

## 2017-05-16 2017-05-17 2017-05-18 2017-05-19 2017-05-20 2017-05-21 
##   334.5000   439.2000   349.2000   345.2000   419.5000   352.9000 
## 2017-05-22 2017-05-23 2017-05-24 2017-05-25 2017-05-26 2017-05-27 
##   372.2000   402.4500   309.5000   382.0000   432.6000   392.6000 
## 2017-05-28 2017-05-29 2017-05-30 2017-05-31 
##   391.0000   405.0000   369.9500   338.2333

length(ag.irregtime)

## [1] 16

####
## Method 2 - date and time component kept
sensor.date1 = strptime(irregular_sensor$X1,
                        '%m/%d/%y %I:%M %p')
sensor.date1

##  [1] "2017-05-16 10:34:00 IST" "2017-05-17 15:23:00 IST"
##  [3] "2017-05-17 20:45:00 IST" "2017-05-18 03:23:00 IST"
##  [5] "2017-05-18 12:34:00 IST" "2017-05-19 11:34:00 IST"
##  [7] "2017-05-20 12:34:00 IST" "2017-05-21 12:34:00 IST"
##  [9] "2017-05-22 17:45:00 IST" "2017-05-22 06:02:00 IST"
## [11] "2017-05-23 04:45:00 IST" "2017-05-23 12:34:00 IST"
## [13] "2017-05-24 02:35:00 IST" "2017-05-25 04:27:00 IST"
## [15] "2017-05-26 15:39:00 IST" "2017-05-27 06:29:00 IST"
## [17] "2017-05-28 07:29:00 IST" "2017-05-29 05:49:00 IST"
## [19] "2017-05-30 07:49:00 IST" "2017-05-30 08:34:00 IST"
## [21] "2017-05-30 13:37:00 IST" "2017-05-30 15:45:00 IST"
## [23] "2017-05-31 05:37:00 IST" "2017-05-31 08:38:00 IST"
## [25] "2017-05-31 16:45:00 IST"

irreg.dates1 = data.frame(1:25) %>% 
  mutate(date = as.Date(sensor.date1),
         measurement = irregular_sensor$X2) %>% 
  select(date, measurement) %>% 
  group_by(date) %>% 
  mutate(measurement = mean(measurement)) %>% 
  arrange(date) %>% 
  distinct()

## Error in select(., date, measurement): unused arguments (date, measurement)

irreg.dates1

## Error in eval(expr, envir, enclos): object 'irreg.dates1' not found

# Creating the zoo object
irreg.dates1 = zoo(irregular_sensor$X2, 
                   order.by = sensor.date1)
irreg.dates1

## 2017-05-16 10:34:00 2017-05-17 15:23:00 2017-05-17 20:45:00 
##               334.5               385.9               492.5 
## 2017-05-18 03:23:00 2017-05-18 12:34:00 2017-05-19 11:34:00 
##               325.8               372.6               345.2 
## 2017-05-20 12:34:00 2017-05-21 12:34:00 2017-05-22 06:02:00 
##               419.5               352.9               392.5 
## 2017-05-22 17:45:00 2017-05-23 04:45:00 2017-05-23 12:34:00 
##               351.9               401.3               403.6 
## 2017-05-24 02:35:00 2017-05-25 04:27:00 2017-05-26 15:39:00 
##               309.5               382.0               432.6 
## 2017-05-27 06:29:00 2017-05-28 07:29:00 2017-05-29 05:49:00 
##               392.6               391.0               405.0 
## 2017-05-30 07:49:00 2017-05-30 08:34:00 2017-05-30 13:37:00 
##               392.5               372.5               312.7 
## 2017-05-30 15:45:00 2017-05-31 05:37:00 2017-05-31 08:38:00 
##               402.1               305.1               402.5 
## 2017-05-31 16:45:00 
##               307.1

plot(irreg.dates1, type = "l")

autoplot(irreg.dates1) +
  xlab("") +
  ylab("Measurement")

# Regularizing with aggregate
ag.irregtime1 = aggregate(irreg.dates1, #data on which to aggregate
                          as.Date, #function on the column by which to aggregate 
                          mean) #aggregating function
ag.irregtime1

## 2017-05-16 2017-05-17 2017-05-18 2017-05-19 2017-05-20 2017-05-21 
##   334.5000   439.2000   349.2000   345.2000   419.5000   352.9000 
## 2017-05-22 2017-05-23 2017-05-24 2017-05-25 2017-05-26 2017-05-27 
##   372.2000   402.4500   309.5000   382.0000   432.6000   392.6000 
## 2017-05-28 2017-05-29 2017-05-30 2017-05-31 
##   391.0000   405.0000   369.9500   338.2333

plot(ag.irregtime1) # plotting the regular zoo object directly

autoplot(ag.irregtime1) +
  xlab("") +
  ylab("Measurement")

myts = ts(ag.irregtime1) # converting to a standard ts, the days start at 1
plot(myts)

autoplot(ag.irregtime1) +
  xlab("") +
  ylab("Measurement")

### Working with Missing Data and Outliers 
## Import ts.NAandOutliers.csv
mydata <- read_csv("./Data/ts-NAandOutliers.csv")

## Warning: Missing column names filled in: 'X1' [1]

## Parsed with column specification:
## cols(
##   X1 = col_double(),
##   mydata = col_double()
## )

# Convert the 2nd column to a simple ts without frequency
myts = ts(mydata$mydata)
myts

## Time Series:
## Start = 1 
## End = 250 
## Frequency = 1 
##   [1]  32.801464  42.465485         NA  32.204058  55.557647  33.050864
##   [7]  43.401620  37.768318  22.844180  36.428877  28.496485  59.037881
##  [13]  36.544163  26.668135  41.325626  28.913199  38.595417  31.341447
##  [19]  34.547023         NA  30.499324  49.391323  43.976004  22.162741
##  [25]  19.439525  41.892407  30.321857  32.899878  17.686235  10.332791
##  [31]  31.612958  40.011275  35.378517  46.167222  26.903207  36.304821
##  [37]  23.408770  42.785841  31.919674  37.571226  33.907485  17.698917
##  [43]  19.931775  23.971169 999.000000  32.853670  33.012320  47.893249
##  [49]  33.961104  40.826518  34.389579  27.210322  41.815827         NA
##  [55]  49.711080  37.246486  34.472507  27.554913  37.976930  24.503481
##  [61]  33.941547  28.582326  17.945402  40.335543  32.103075  15.609346
##  [67]  38.637130  58.877558  42.178769  34.075469  29.208206  20.409934
##  [73]  23.682860  49.014566  59.160903  24.994359  37.321672  11.830421
##  [79]  49.907975  33.288427  25.900307  34.661099  38.170951  30.246685
##  [85]  45.001326  36.082827  38.969588  24.260726   8.619401  33.933167
##  [91]  30.158056  32.211135  46.688584  36.399098  27.266510  39.706101
##  [97]  48.560701 999.000000  31.011612  33.565184  41.850476  45.780926
## [103]  21.679404  32.340497  55.904896  17.349895  32.994516  36.155426
## [109]  47.089342  33.955275  36.563838  18.773382  28.077605  40.483324
## [115]  41.341771  32.907839  59.604911  20.989279  46.886734  53.931163
## [121]  44.662468  43.125045  25.800244  22.833920  51.397357  34.775922
## [127]  50.922532  36.430258  32.975690  37.659017  48.006323  49.901919
## [133]  20.619643  24.895206   4.682543  17.049461  45.618543  28.288209
## [139]  50.446258  34.983971  38.847283  32.301493  46.044574  21.739473
## [145]  16.457915  36.157602  35.773314  23.368300  34.220736  39.443674
## [151]  26.074044  28.599269  45.410516         NA  30.585004  23.405284
## [157]  36.949438  17.647508  22.991044  40.388899  30.654440  49.261182
## [163]  31.215505  30.462442  41.294423  28.046393  24.925970  23.934094
## [169]  41.690112  46.226476  40.741721 999.000000  26.455048  12.766929
## [175]  38.133315  61.653241  11.474054  41.835335  34.572898  59.921615
## [181]  53.072688  51.889866  43.700888  33.639492  24.988901  27.019376
## [187]  12.774882  21.001551  18.133113  48.563807  64.633075  29.362668
## [193]  45.478245  34.152629  50.573869  43.564419  57.644084  20.785408
## [199]  39.811707  51.854335  33.351763  23.242107  34.351879  28.113792
## [205]  33.952911  35.372668  38.059841  40.818440  45.768335  39.272128
## [211] 999.000000  36.665537  50.282893  34.561040  46.040830  39.936308
## [217]  39.873144  51.126396   2.683472  51.667975  41.336229  43.090450
## [223]  24.686842  52.300908  37.379943  19.043254  43.512121  54.236360
## [229]  33.557160  33.851597         NA  36.149460  32.985037  31.422766
## [235]  36.574851  29.648483  18.290954  38.154075  34.446452  12.037743
## [241]  23.581530  36.968395  50.747174  37.981389  28.693203  32.396009
## [247]  41.325484  30.017571  14.818111  45.403854

# Checking for NAs and outliers
summary(myts)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   2.683  28.078  34.573  50.710  42.465 999.000       5

plot(myts)

autoplot(myts) +
  xlab("") +
  ylab("")

# Automatic detection of outliers
myts1 = tsoutliers(myts)
myts1

## $index
## [1]  45  98 172 211
## 
## $replacements
## [1] 28.41242 39.78616 33.59838 37.96883

plot(myts)

autoplot(myts) +
  xlab("") +
  ylab("")

# Missing data handling with zoo
myts.NAlocf = na.locf(myts)
myts.NAfill = na.fill(myts, 33)

# Tip: na.trim to get rid of NAs at the beginning or end of dataset
# Standard NA method in package forecast
myts.NAinterp = na.interp(myts)

# Cleaning NA and outliers with forecast package
mytsclean = tsclean(myts)
plot(mytsclean)

autoplot(mytsclean) +
  xlab("") +
  ylab("")

summary(mytsclean)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.683  28.157  34.567  35.025  41.830  64.633