title : “Time Serries Forecasting”

author : “Recrafted by bambangpe”

output : html_document

Upload Library

Data Preparation

monthly_milk <- read.csv("D:/data/monthly_milk.csv")
# Milk production per cow per month
head(monthly_milk)

##        month milk_prod_per_cow_kg
## 1 1962-01-01               265.05
## 2 1962-02-01               252.45
## 3 1962-03-01               288.00
## 4 1962-04-01               295.20
## 5 1962-05-01               327.15
## 6 1962-06-01               313.65

daily_milk <- read.csv("D:/data/daily_milk.csv")
# Milk production per cow per milking
head(daily_milk)

##             date_time milk_prod_per_cow_kg
## 1 1975-01-01 05:00:00             11.21745
## 2 1975-01-01 17:00:00             10.67182
## 3 1975-01-02 05:00:00             10.90791
## 4 1975-01-02 17:00:00             11.03970
## 5 1975-01-03 05:00:00             12.53303
## 6 1975-01-03 17:00:00             10.69446

class(monthly_milk)

## [1] "data.frame"

Coerce monthly_milk to Date class

class(monthly_milk$month)

## [1] "factor"

monthly_milk$month_date <- as.Date(monthly_milk$month, format = "%Y-%m-%d")
# Check it worked
class(monthly_milk$month_date) 

## [1] "Date"

Experiment with format()

format(monthly_milk$month_date, format = "%Y-%B-%u")
format(monthly_milk$month_date, format = "%Y-%B-%u")

class(format(monthly_milk$month_date, format = "%Y-%B-%u"))  # class is no longer `Date`

## [1] "character"

Coerce daily milk to POSIXct class

head(daily_milk)

##             date_time milk_prod_per_cow_kg
## 1 1975-01-01 05:00:00             11.21745
## 2 1975-01-01 17:00:00             10.67182
## 3 1975-01-02 05:00:00             10.90791
## 4 1975-01-02 17:00:00             11.03970
## 5 1975-01-03 05:00:00             12.53303
## 6 1975-01-03 17:00:00             10.69446

head(daily_milk)

##             date_time milk_prod_per_cow_kg
## 1 1975-01-01 05:00:00             11.21745
## 2 1975-01-01 17:00:00             10.67182
## 3 1975-01-02 05:00:00             10.90791
## 4 1975-01-02 17:00:00             11.03970
## 5 1975-01-03 05:00:00             12.53303
## 6 1975-01-03 17:00:00             10.69446

class(daily_milk$date_time)

## [1] "factor"

daily_milk$date_time_posix <- as.POSIXct(daily_milk$date_time, format = "%Y-%m-%d %H:%M:%S")

Check it worked

class(daily_milk$date_time_posix)

## [1] "POSIXct" "POSIXt"

Create badly formatted dates

## [1] "01/Jan/1962-1" "01/Feb/1962-4" "01/Mar/1962-4" "01/Apr/1962-7"
## [5] "01/May/1962-2" "01/Jun/1962-5"

## [1] "character"

Coerce bad date to good date

monthly_milk$good_date <- as.Date(monthly_milk$bad_date, format = "%d/%b/%Y-%u")
head(monthly_milk$good_date)

## [1] "1962-01-01" "1962-02-01" "1962-03-01" "1962-04-01" "1962-05-01"
## [6] "1962-06-01"

class(monthly_milk$good_date)

## [1] "Date"

Plot time series data

Using scale_x_datetime

Viewing trend using loess smooth

## `geom_smooth()` using formula 'y ~ x'

Explore seasonal trends

Extract month and year and store in new column

monthly_milk$year <- format(monthly_milk$month_date, format = "%Y")
monthly_milk$month_num <- format(monthly_milk$month_date, format = "%m")

Create a colour palette using the colortools package

Plot months

Using ts objects to decompose trends

Transform to ts class

monthly_milk_ts <- ts(monthly_milk$milk_prod, start = 1962, end = 1975, freq = 12)  
# Specify start and end year, measurement frequency (monthly = 12)

Decompose using stl()

monthly_milk_stl <- stl(monthly_milk_ts, s.window = "period")

Generate plots

monthplot(monthly_milk_ts, choice = "seasonal")  # variation in milk production for each month

seasonplot(monthly_milk_ts) #ok

Forecasting

# Split data into testing and training
monthly_milk_model <- window(x = monthly_milk_ts, start=c(1962), end=c(1970))
monthly_milk_test <- window(x = monthly_milk_ts, start=c(1970))

Creating model objects of each type of ets model

milk_ets_auto <- ets(monthly_milk_model)
milk_ets_mmm <- ets(monthly_milk_model, model = "MMM")
milk_ets_zzz<- ets(monthly_milk_model, model = "ZZZ")
milk_ets_mmm_damped <- ets(monthly_milk_model, model = "MMM", damped = TRUE)

Creating forecast objects

milk_ets_fc <- forecast(milk_ets_auto, h = 60)  # `h = 60` means that the forecast will be 60 time periods long, in our case a time period is one month
milk_ets_mmm_fc <- forecast(milk_ets_mmm, h = 60)
milk_ets_zzz_fc <- forecast(milk_ets_zzz, h = 60)
milk_ets_mmm_damped_fc <- forecast(milk_ets_mmm_damped, h = 60)

Convert forecasts to data frames

Combining into one data frame

forecast_all <- rbind(milk_ets_fc_df, milk_ets_mmm_fc_df, milk_ets_zzz_fc_df, milk_ets_mmm_damped_fc_df)

Plotting with ggplot

Comparing accuracy of forecasts

accuracy(milk_ets_fc, monthly_milk_test) accuracy(milk_ets_mmm_fc, monthly_milk_test) accuracy(milk_ets_zzz_fc, monthly_milk_test) accuracy(milk_ets_mmm_damped_fc, monthly_milk_test)

Extracting forecast estimates

milk_ets_fc_df %>%
    filter(Month == "Jan 1975") %>%
    select(Month, Point_Forecast)

##      Month Point_Forecast
## 1 Jan 1975       383.7897

milk_ets_zzz_fc_df %>%
    filter(Month == "Jan 1975") %>%
    select(Month, Point_Forecast)

##      Month Point_Forecast
## 1 Jan 1975       383.7897