Baltmore city annual water use, liters per capita per day, 1885-1968

water = read.csv("baltmore-city-annual-water-use-l.csv")
library(fpp)

## Loading required package: forecast

## Warning: package 'forecast' was built under R version 3.3.2

## Loading required package: fma

## Warning: package 'fma' was built under R version 3.3.2

## Loading required package: expsmooth

## Loading required package: lmtest

## Warning: package 'lmtest' was built under R version 3.3.2

## Loading required package: zoo

## Warning: package 'zoo' was built under R version 3.3.2

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## Loading required package: tseries

## Warning: package 'tseries' was built under R version 3.3.2

library(forecast)
library(xts)

## Warning: package 'xts' was built under R version 3.3.2

library(tseries)

str(water)

## 'data.frame':    79 obs. of  2 variables:
##  $ Year                                                                : int  1885 1886 1887 1888 1889 1890 1891 1892 1893 1894 ...
##  $ Baltmore.city.annual.water.use..liters.per.capita.per.day..1885.1968: int  356 386 397 397 413 458 485 344 390 360 ...

head(water)

##   Year
## 1 1885
## 2 1886
## 3 1887
## 4 1888
## 5 1889
## 6 1890
##   Baltmore.city.annual.water.use..liters.per.capita.per.day..1885.1968
## 1                                                                  356
## 2                                                                  386
## 3                                                                  397
## 4                                                                  397
## 5                                                                  413
## 6                                                                  458

tail (water)

##    Year
## 74 1958
## 75 1959
## 76 1960
## 77 1961
## 78 1962
## 79 1963
##    Baltmore.city.annual.water.use..liters.per.capita.per.day..1885.1968
## 74                                                                  602
## 75                                                                  594
## 76                                                                  587
## 77                                                                  587
## 78                                                                  625
## 79                                                                  613

names(water)[2]<-c("use")

train = water[1:63,]
test = water[64:79,]

time series

ts.water = ts(water$use, frequency = 1, start = c(1885,1))
plot(ts.water)

# 63 years
ts.train =ts(train$use,frequency=1,start=c(1885,1))
plot(ts.train)

# 16 years
ts.test = ts(test$use, frequency=1, start=c(1948))
plot(ts.test)

Models: the data are not seasonal so, no seasonal models were used

fit1 = ses(ts.train)
fit2 = holt(ts.train)
fit3 = holt(ts.train,exponential=TRUE)
fit4 = holt(ts.train,damped=TRUE)
fit5 = holt(ts.train,exponential=TRUE,damped=TRUE)

# Results for first model:
fit1$model

## Simple exponential smoothing 
## 
## Call:
##  ses(y = ts.train) 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
## 
##   Initial states:
##     l = 355.9419 
## 
##   sigma:  35.6975
## 
##      AIC     AICc      BIC 
## 717.4776 717.8844 723.9070

accuracy(fit1) # training set

##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 4.144218 35.69748 23.79471 0.5850507 4.952862 0.9841707
##                    ACF1
## Training set 0.01900003

accuracy(fit1,ts.test) # test set

##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set   4.144218 35.69748 23.79471  0.5850507 4.952862 0.9841707
## Test set     -42.401501 47.24963 42.40150 -7.5208965 7.520897 1.7537646
##                    ACF1 Theil's U
## Training set 0.01900003        NA
## Test set     0.55384686  3.160426

ME1=  -42.401501
RMSE1 = 47.24963    
MAE1 = 42.40150   
MPE1= -7.5208965
MAPE1= 7.520897
MASE1= 1.7537646   
     
# Results for second model:
fit2$model

## Holt's method 
## 
## Call:
##  holt(y = ts.train) 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
##     beta  = 1e-04 
## 
##   Initial states:
##     l = 396.5706 
##     b = 3.8617 
## 
##   sigma:  35.8949
## 
##      AIC     AICc      BIC 
## 722.1726 723.2253 732.8883

accuracy(fit2) # training set

##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -0.3601994 35.89495 23.88767 -0.3949746 5.033001 0.9880158
##                     ACF1
## Training set 0.005778681

accuracy(fit2,ts.test) # test set

##                       ME     RMSE      MAE         MPE      MAPE      MASE
## Training set  -0.3601994 35.89495 23.88767  -0.3949746  5.033001 0.9880158
## Test set     -63.6286874 68.90130 63.62869 -11.2399795 11.239979 2.6317402
##                     ACF1 Theil's U
## Training set 0.005778681        NA
## Test set     0.590756262   4.59837

ME2=   -63.6286874
RMSE2 = 68.90130    
MAE2 = 63.62869   
MPE2= -11.2399795
MAPE2= 11.239979
MASE2= 2.6317402

      
# Results for third model:
fit3$model

## Holt's method with exponential trend 
## 
## Call:
##  holt(y = ts.train, exponential = TRUE) 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
##     beta  = 1e-04 
## 
##   Initial states:
##     l = 397.21 
##     b = 1.0047 
## 
##   sigma:  0.0738
## 
##      AIC     AICc      BIC 
## 722.3212 723.3738 733.0368

accuracy(fit3) # training set

##                    ME    RMSE      MAE         MPE     MAPE     MASE
## Training set 1.199121 36.0041 24.20351 -0.06419636 5.093922 1.001079
##                     ACF1
## Training set 0.005386954

accuracy(fit3,ts.test) # test set

##                      ME     RMSE      MAE          MPE      MAPE     MASE
## Training set   1.199121 36.00410 24.20351  -0.06419636  5.093922 1.001079
## Test set     -58.527237 63.55173 58.52724 -10.34605296 10.346053 2.420740
##                     ACF1 Theil's U
## Training set 0.005386954        NA
## Test set     0.580782417  4.243352

ME3=  -58.527237
RMSE3 = 63.55173    
MAE3 = 58.52724   
MPE3= -10.34605296
MAPE3= 10.346053
MASE3= 2.420740
      
            
# Results for fourth model:
fit4$model

## Damped Holt's method 
## 
## Call:
##  holt(y = ts.train, damped = TRUE) 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
##     beta  = 1e-04 
##     phi   = 0.98 
## 
##   Initial states:
##     l = 395.1801 
##     b = 5.2393 
## 
##   sigma:  35.903
## 
##      AIC     AICc      BIC 
## 724.2010 725.7010 737.0598

accuracy(fit4) # training set

##                     ME     RMSE    MAE        MPE     MAPE      MASE
## Training set 0.5882346 35.90303 23.927 -0.2209718 5.033189 0.9896423
##                     ACF1
## Training set 0.007069687

accuracy(fit4,ts.test) # test set

##                       ME     RMSE    MAE        MPE     MAPE      MASE
## Training set   0.5882346 35.90303 23.927 -0.2209718 5.033189 0.9896423
## Test set     -49.8780045 54.66056 49.878 -8.8310543 8.831054 2.0629995
##                     ACF1 Theil's U
## Training set 0.007069687        NA
## Test set     0.563391484  3.652957

ME4=  -49.8780045
RMSE4 = 54.66056    
MAE4 = 49.878   
MPE4= -8.8310543
MAPE4= 8.831054
MASE4= 2.0629995

# Results for fifth model:
fit5$model

## Damped Holt's method with exponential trend 
## 
## Call:
##  holt(y = ts.train, damped = TRUE, exponential = TRUE) 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
##     beta  = 1e-04 
##     phi   = 0.8653 
## 
##   Initial states:
##     l = 395.0726 
##     b = 1.0013 
## 
##   sigma:  0.0744
## 
##      AIC     AICc      BIC 
## 724.8878 726.3878 737.7466

accuracy(fit5) # training set

##                    ME     RMSE     MAE       MPE     MAPE     MASE
## Training set 3.467551 36.03582 24.4017 0.3971757 5.124481 1.009276
##                     ACF1
## Training set 0.006067191

accuracy(fit5,ts.test) # test set

##                      ME     RMSE      MAE        MPE     MAPE     MASE
## Training set   3.467551 36.03582 24.40170  0.3971757 5.124481 1.009276
## Test set     -42.420324 47.26807 42.42032 -7.5241979 7.524198 1.754543
##                     ACF1 Theil's U
## Training set 0.006067191        NA
## Test set     0.553860084  3.161652

ME5=  -42.420324
RMSE5 = 47.26807    
MAE5 = 42.42032   
MPE5= -7.5241979
MAPE5= 7.524198
MASE5= 1.754543

##Comparing errors
error.matrix=matrix(c(ME1, RMSE1,MAE1,MPE1, MAPE1, MASE1,ME2, RMSE2,MAE2,MPE2, MAPE2, MASE2,ME3, RMSE3,MAE3,MPE3, MAPE3, MASE3,ME4, RMSE4,MAE4,MPE4, MAPE4, MASE4,ME5, RMSE5,MAE5,MPE5, MAPE5, MASE5),ncol=5)
colnames(error.matrix)=c("SES","Holt's","Exponential", "Additive", "Multiplicative")
rownames(error.matrix)=c("ME","RMSE","MAE","MPE","MAPE", "MASE")
error.table<-as.table(error.matrix)
error.table

##             SES     Holt's Exponential   Additive Multiplicative
## ME   -42.401501 -63.628687  -58.527237 -49.878005     -42.420324
## RMSE  47.249630  68.901300   63.551730  54.660560      47.268070
## MAE   42.401500  63.628690   58.527240  49.878000      42.420320
## MPE   -7.520897 -11.239980  -10.346053  -8.831054      -7.524198
## MAPE   7.520897  11.239979   10.346053   8.831054       7.524198
## MASE   1.753765   2.631740    2.420740   2.063000       1.754543

Plot comparing forecasting performance of non-seasonal methods

plot(fit3, type="o", ylab="Annual water use (liters per capita per day)", 
  flwd=1)
lines(window(ts.train,start=1885),type="o")
lines(fit1$mean,col=2)
lines(fit2$mean,col=3)
lines(fit4$mean,col=5)
lines(fit5$mean,col=6)
legend("topleft", lty=1, pch=1, col=1:6,
    c("Data","SES","Holt's","Exponential",
      "Additive Damped","Multiplicative Damped"))

Discussion 3

Baltmore city annual water use, liters per capita per day, 1885-1968

time series

Models: the data are not seasonal so, no seasonal models were used

Plot comparing forecasting performance of non-seasonal methods