R Notebook

1.a) Visualizing the whole data,dividing the time series data into different components and plotting the different components individually:

library(forecast)
library(tseries)

Data Import

setwd("C:/Users/Divya Vignesh/Desktop/BABI/Time Series")

The working directory was changed to C:/Users/Divya Vignesh/Desktop/BABI/Time Series inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

data <- read.csv("fancy.csv",header = TRUE)

Summary

summary(data)

 SalesSouvenir   
 Min.   :  1665  
 1st Qu.:  5884  
 Median :  8772  
 Mean   : 14316  
 3rd Qu.: 16889  
 Max.   :104661  

str(data)

'data.frame':   84 obs. of  1 variable:
 $ SalesSouvenir: num  1665 2398 2841 3547 3753 ...

View(data)
head(data)

Framing Time series data

?ts
tsdata <- ts(data=data,start=c(1987,1),end=c(1993,12),frequency=12)
tsdata

           Jan       Feb       Mar       Apr       May       Jun       Jul
1987   1664.81   2397.53   2840.71   3547.29   3752.96   3714.74   4349.61
1988   2499.81   5198.24   7225.14   4806.03   5900.88   4951.34   6179.12
1989   4717.02   5702.63   9957.58   5304.78   6492.43   6630.80   7349.62
1990   5921.10   5814.58  12421.25   6369.77   7609.12   7224.75   8121.22
1991   4826.64   6470.23   9638.77   8821.17   8722.37  10209.48  11276.55
1992   7615.03   9849.69  14558.40  11587.33   9332.56  13082.09  16732.78
1993  10243.24  11266.88  21826.84  17357.33  15997.79  18601.53  26155.15
           Aug       Sep       Oct       Nov       Dec
1987   3566.34   5021.82   6423.48   7600.60  19756.21
1988   4752.15   5496.43   5835.10  12600.08  28541.72
1989   8176.62   8573.17   9690.50  15151.84  34061.01
1990   7979.25   8093.06   8476.70  17914.66  30114.41
1991  12552.22  11637.39  13606.89  21822.11  45060.69
1992  19888.61  23933.38  25391.35  36024.80  80721.71
1993  28586.52  30505.41  30821.33  46634.38 104660.67

Plotting the framed Time series data

# the stationary signal and ACF
plot(tsdata,bty="l",
     type='l',
     xlab = "Year",
     ylab = "Sales of sovenier data",
     main = "Fancy sovenier Sales")
sm <- ma(tsdata, order=12) # 12 month moving average
lines(sm, col="red")

Inferences :

1.We have sales for all the months from the year from 1987 to 1993

2.As per the plots and the Data Framed we could see that there is gradual increase in the sales as the years go by

3.There seems to be a spike increase in sales in the month of December every year which seems to be a peak season for sales

4.And the small mountain peaks where in the month of march every year where the sales slightly goes up and comes back to normal as per the data

5.We could observe seasonality in which November and December have very high sales

6. There has been slight dip in the year 1991 which can be observed.

par(oma=c(0,0,0,2))

xrange <- c(1987,1993)
yrange <- range(tsdata)
plot(xrange, yrange, type="n", xlab="Date", ylab="Sales",  bty="l", las=1)

colors <- rainbow(12)
linetype <- c(1:12)
plotchar <- c(1:12)

axis(1, at=seq(1995,2001,1), labels=format(seq(1995,2001,1)))

for (i in 1:12) {
  currentMonth <- subset(tsdata, cycle(tsdata)==i)
  lines(seq(1987, 1987+length(currentMonth)-1,1), currentMonth, type="b", lwd=1,
      lty=linetype[i], col=colors[i], pch=plotchar[i])
}

title("Souvenir Sales Broken Out by Month")

legend(2001.35,100000, 1:12, cex=0.8, col=colors, pch=plotchar, lty=linetype, title="Month", xpd=NA)

require(graphics)

plot(stlper<-stl(tsdata[,1], "per"))

summary(stlper)

 Call:
 stl(x = tsdata[, 1], s.window = "per")

 Time.series components:
    seasonal            trend            remainder         
 Min.   :-6875.91   Min.   : 5835.38   Min.   :-18530.127  
 1st Qu.:-5061.09   1st Qu.: 8151.91   1st Qu.: -1843.078  
 Median :-2911.02   Median :10639.86   Median :   400.246  
 Mean   :    0.00   Mean   :14540.62   Mean   :  -225.036  
 3rd Qu.:-1500.37   3rd Qu.:16922.29   3rd Qu.:  1887.459  
 Max.   :32250.85   Max.   :45025.84   Max.   : 27383.973  
 IQR:
     STL.seasonal STL.trend STL.remainder data 
      3561         8770      3731         11004
   %  32.4         79.7      33.9         100.0

 Weights: all == 1

 Other components: List of 5
 $ win  : Named num [1:3] 841 19 13
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 85 2 2
 $ inner: int 2
 $ outer: int 0

plot(stlswin<-stl(tsdata[,1], s.window = 7, t.window = 50, t.jump = 1))

summary(stlswin)

 Call:
 stl(x = tsdata[, 1], s.window = 7, t.window = 50, t.jump = 1)

 Time.series components:
    seasonal            trend            remainder        
 Min.   :-9161.09   Min.   : 4856.96   Min.   :-9305.964  
 1st Qu.:-3798.68   1st Qu.: 8144.81   1st Qu.:-1689.841  
 Median :-2440.99   Median :11065.64   Median : -223.507  
 Mean   :  231.84   Mean   :14447.32   Mean   : -363.577  
 3rd Qu.:-1220.21   3rd Qu.:19184.16   3rd Qu.:  875.385  
 Max.   :43249.71   Max.   :35695.01   Max.   :25715.952  
 IQR:
     STL.seasonal STL.trend STL.remainder data 
      2578        11039      2565         11004
   %  23.4        100.3      23.3         100.0

 Weights: all == 1

 Other components: List of 5
 $ win  : Named num [1:3] 7 50 13
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 1 1 2
 $ inner: int 2
 $ outer: int 0

plot(stllc <- stl(log(tsdata[,1]), s.window = 21))

summary(stllc)

 Call:
 stl(x = log(tsdata[, 1]), s.window = 21)

 Time.series components:
    seasonal              trend             remainder          
 Min.   :-0.6435919   Min.   : 8.108900   Min.   :-0.29879751  
 1st Qu.:-0.2681912   1st Qu.: 8.838632   1st Qu.:-0.06822801  
 Median :-0.0653205   Median : 9.128885   Median : 0.00090540  
 Mean   : 0.0002944   Mean   : 9.218226   Mean   : 0.00147858  
 3rd Qu.: 0.0521100   3rd Qu.: 9.625153   3rd Qu.: 0.08686454  
 Max.   : 1.2852383   Max.   :10.304718   Max.   : 0.28224802  
 IQR:
     STL.seasonal STL.trend STL.remainder data  
     0.3203       0.7865    0.1551        1.0542
   %  30.4         74.6      14.7         100.0 

 Weights: all == 1

 Other components: List of 5
 $ win  : Named num [1:3] 21 21 13
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 3 3 2
 $ inner: int 2
 $ outer: int 0

## linear trend, strict period.
plot(stllcwin<-stl(log(tsdata[,1]), s.window = "per", t.window = 1000))
summary(stllcwin)

 Call:
 stl(x = log(tsdata[, 1]), s.window = "per", t.window = 1000)

 Time.series components:
    seasonal              trend             remainder         
 Min.   :-0.6649476   Min.   : 8.291322   Min.   :-0.4154458  
 1st Qu.:-0.2615064   1st Qu.: 8.755716   1st Qu.:-0.1247289  
 Median :-0.0651962   Median : 9.220109   Median : 0.0059606  
 Mean   : 0.0000000   Mean   : 9.220109   Mean   :-0.0001103  
 3rd Qu.: 0.0441778   3rd Qu.: 9.684502   3rd Qu.: 0.1123122  
 Max.   : 1.2982109   Max.   :10.148896   Max.   : 0.3872554  
 IQR:
     STL.seasonal STL.trend STL.remainder data  
     0.3057       0.9288    0.2370        1.0542
   %  29.0         88.1      22.5         100.0 

 Weights: all == 1

 Other components: List of 5
 $ win  : Named num [1:3] 841 1000 13
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 85 100 2
 $ inner: int 2
 $ outer: int 0

## Two STL plotted side by side :
stmd <- stl(tsdata[,1], s.window = "per") # non-robust
summary(stmd)

 Call:
 stl(x = tsdata[, 1], s.window = "per")

 Time.series components:
    seasonal            trend            remainder         
 Min.   :-6875.91   Min.   : 5835.38   Min.   :-18530.127  
 1st Qu.:-5061.09   1st Qu.: 8151.91   1st Qu.: -1843.078  
 Median :-2911.02   Median :10639.86   Median :   400.246  
 Mean   :    0.00   Mean   :14540.62   Mean   :  -225.036  
 3rd Qu.:-1500.37   3rd Qu.:16922.29   3rd Qu.:  1887.459  
 Max.   :32250.85   Max.   :45025.84   Max.   : 27383.973  
 IQR:
     STL.seasonal STL.trend STL.remainder data 
      3561         8770      3731         11004
   %  32.4         79.7      33.9         100.0

 Weights: all == 1

 Other components: List of 5
 $ win  : Named num [1:3] 841 19 13
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 85 2 2
 $ inner: int 2
 $ outer: int 0

summary(stmR <- stl(tsdata[,1], s.window = "per", robust = TRUE))

 Call:
 stl(x = tsdata[, 1], s.window = "per", robust = TRUE)

 Time.series components:
    seasonal             trend            remainder        
 Min.   :-5837.471   Min.   : 6746.77   Min.   :-18332.91  
 1st Qu.:-4064.275   1st Qu.: 9030.87   1st Qu.: -1161.33  
 Median :-2796.039   Median :11103.34   Median :   -52.07  
 Mean   :    0.000   Mean   :14314.28   Mean   :     1.31  
 3rd Qu.:-2090.811   3rd Qu.:16831.90   3rd Qu.:   650.98  
 Max.   :29650.647   Max.   :36130.75   Max.   : 38879.27  
 IQR:
     STL.seasonal STL.trend STL.remainder data 
      1973         7801      1812         11004
   %  17.9         70.9      16.5         100.0

 Weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.6752  0.9452  0.7342  0.9879  1.0000 

 Other components: List of 5
 $ win  : Named num [1:3] 841 19 13
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 85 2 2
 $ inner: int 1
 $ outer: int 15

op <- par(mar = c(0, 4, 0, 3), oma = c(5, 0, 4, 0), mfcol = c(4, 2))

plot(stmd, set.pars = NULL, labels  =  NULL,
     main = "stl(my_data.ts[,1], s.w = \"per\",  robust = FALSE / TRUE )")
plot(stmR, set.pars = NULL)

stmR $ weights

 [1] 0.95320350 0.99124572 0.67983843 0.94129633 0.99956318 0.99618458
 [7] 0.99278201 0.80519116 0.99505138 0.98778808 0.00000000 0.00000000
[13] 0.99441119 0.86929698 0.98097476 0.98733036 0.96649001 0.99945635
[19] 0.98830206 0.82738770 0.94902843 0.88984748 0.71358774 0.00000000
[25] 0.90682184 0.94839133 0.66123688 0.98393780 0.99030618 0.99999516
[31] 0.98762003 0.99987127 0.98653438 0.94007156 0.89174838 0.00000000
[37] 0.96116480 0.98236562 0.25947144 0.99366711 0.99880493 0.99968030
[43] 0.99986773 0.97590523 0.98561888 0.95535944 0.88710044 0.00000000
[49] 0.89000369 0.93969350 0.96117156 0.98166800 0.94369446 0.96949629
[55] 0.99243104 0.93186964 0.99274004 0.94671769 0.85584518 0.99967776
[61] 0.59141497 0.79960726 0.99995806 0.87683494 0.00000000 0.77270481
[67] 0.99950586 0.60577289 0.00000000 0.00000000 0.00000000 0.00000000
[73] 0.00000000 0.00000000 0.90019832 0.05403848 0.00000000 0.00000000
[79] 0.98114218 0.99044908 0.91665144 0.97192358 0.00000000 0.00000000

# mark the 'outliers' :
(iO <- which(stmR $ weights  < 1e-8)) # 10 were considered outliers

 [1] 11 12 24 36 48 65 69 70 71 72 73 74 77 78 83 84

sts <- stmR$time.series
points(time(sts)[iO], 0.8* sts[,"remainder"][iO], pch = 4, col = "red")
par(op)   # reset

Decomposing Time Series

decompose_multi<-decompose(tsdata,type='m')
plot(decompose(tsdata,type='m'))

plot(decompose_multi)

?decompose
decomposets<- decompose(tsdata)
plot(decomposets)

Seasonal Plot

seasonplot(tsdata)

Here with respect to season same thing as infered before that march and Decemeber has some high sales compared to other months.

On comparing multiplicative and additive time series,additive time series we are able to observe good pattern which can be useful for future time series prediction. Also we will transform them into log scale and continue our time series forecasting.

Transforming time series data

Taking Log on the time series for using additive model.Additive model is used when the variance of the time series doesn’t change over different values of the time series.

logtsdata<-log(tsdata)
plot(logtsdata)

Decomposing Time Series

logdecomposets<- decompose(logtsdata,type = c("additive", "multiplicative"), filter = NULL)
plot(logdecomposets)

Seasonal a. March & December are having higher sales with December being highest in sales. Trend a. From the graph we are able to observe a positive trend with a spike after 1991. Random a. From the graph we are able to observe a positive trend with a spike after 1991.

Seasonal Plot on transformed data

seasonplot(logtsdata)

Here with respect to season same thing as infered before that march and Decemeber has some high sales compared to other months.

1.b)Time series stationarity check - KPSS Test

kpss.test(logtsdata)

p-value smaller than printed p-value

    KPSS Test for Level Stationarity

data:  logtsdata
KPSS Level = 1.7909, Truncation lag parameter = 3, p-value = 0.01

Observation:

We could see that the p-value of kpss is 0.01 which tells us that the series is non stationary

2.Splitting data into Test/Train or Dev/holdout

logtsdatatrain<-window(logtsdata,start=c(1987,1),end=c(1991,12))

logtsdatatrain

           Jan       Feb       Mar       Apr       May       Jun       Jul
1987  7.417466  7.782194  7.951809  8.173939  8.230300  8.220064  8.377841
1988  7.823970  8.556075  8.885322  8.477627  8.682857  8.507414  8.728931
1989  8.458933  8.648683  9.206089  8.576364  8.778392  8.799481  8.902404
1990  8.686278  8.668124  9.427164  8.759319  8.937103  8.885268  9.002236
1991  8.481906  8.774967  9.173549  9.084910  9.073646  9.231072  9.330481
           Aug       Sep       Oct       Nov       Dec
1987  8.179295  8.521548  8.767715  8.935982  9.891223
1988  8.466352  8.611854  8.671647  9.441458 10.259122
1989  9.009034  9.056393  9.178901  9.625877 10.435909
1990  8.984600  8.998762  9.045077  9.793375 10.312759
1991  9.437653  9.361978  9.518332  9.990679 10.715766

logtsdatatest<-window(logtsdata,start=c(1992,1),end=c(1993,12))

logtsdatatest

           Jan       Feb       Mar       Apr       May       Jun       Jul
1992  8.937879  9.195195  9.585923  9.357668  9.141265  9.478999  9.725125
1993  9.234373  9.329623  9.990896  9.761770  9.680206  9.830999 10.171801
           Aug       Sep       Oct       Nov       Dec
1992  9.897902 10.083029 10.142164 10.491963 11.298763
1993 10.260691 10.325659 10.335962 10.750093 11.558479

tsdataholt <- HoltWinters(logtsdatatrain)

optimization difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH

tsdataholt

Holt-Winters exponential smoothing with trend and additive seasonal component.

Call:
HoltWinters(x = logtsdatatrain)

Smoothing parameters:
 alpha: 0.3629272
 beta : 0
 gamma: 0.8610514

Coefficients:
           [,1]
a    9.57155473
b    0.02996319
s1  -0.64883577
s2  -0.42100919
s3   0.08160850
s4  -0.21696070
s5  -0.21714670
s6  -0.16491424
s7  -0.10613895
s8  -0.09228018
s9  -0.17305267
s10 -0.06456957
s11  0.45893851
s12  1.14554103

plot(tsdataholt)

Observation :

We see from the plot that the Holt-Winters method is predicting right with respect to seasonalpeaks, which occur roughly in December every year

Forecast using Holts winter model

tsdataholtforecast <- forecast(tsdataholt,60)
plot(tsdataholtforecast)

model<-stlf(logtsdatatrain,60)
plot(model)

Observations

The blue line is showing forecast value forecasted by Holt’s Winter model (training), and the Dark and ligth grey areas show 80% and 95% prediction intervals,

Comparison b/n the predicted values and Test data

tsdataholtforecast

logtsdatatest

           Jan       Feb       Mar       Apr       May       Jun       Jul
1992  8.937879  9.195195  9.585923  9.357668  9.141265  9.478999  9.725125
1993  9.234373  9.329623  9.990896  9.761770  9.680206  9.830999 10.171801
           Aug       Sep       Oct       Nov       Dec
1992  9.897902 10.083029 10.142164 10.491963 11.298763
1993 10.260691 10.325659 10.335962 10.750093 11.558479

Testing for accuracy for MAPE

accuracy(tsdataholtforecast,logtsdatatest)

                      ME      RMSE       MAE        MPE     MAPE      MASE
Training set -0.01590626 0.1681845 0.1313730 -0.1789673 1.461293 0.4666977
Test set      0.02907493 0.1916927 0.1603372  0.2205072 1.609411 0.5695922
                   ACF1 Theil's U
Training set -0.1149566        NA
Test set      0.6598762 0.3849737

Inferences

MAPE is 1.46 and 1.60 for training and test data set respectively, This means the model is stable with an accuracy of 1.46% and 1.60% for training and test data set respectively

Predicting for 5 years using the compelete data set

logtsdataholt <- HoltWinters(logtsdata)
logtsdataholt

Holt-Winters exponential smoothing with trend and additive seasonal component.

Call:
HoltWinters(x = logtsdata)

Smoothing parameters:
 alpha: 0.413418
 beta : 0
 gamma: 0.9561275

Coefficients:
           [,1]
a   10.37661961
b    0.02996319
s1  -0.80952063
s2  -0.60576477
s3   0.01103238
s4  -0.24160551
s5  -0.35933517
s6  -0.18076683
s7   0.07788605
s8   0.10147055
s9   0.09649353
s10  0.05197826
s11  0.41793637
s12  1.18088423

plot(logtsdataholt)

tsdataholtfullforecast <- forecast(logtsdataholt,60)
plot(tsdataholtfullforecast)

Checking accuracy

accuracy(tsdataholtfullforecast)

                       ME      RMSE       MAE         MPE     MAPE      MASE
Training set -0.006421623 0.1671448 0.1278046 -0.08896552 1.385112 0.4077539
                   ACF1
Training set 0.04525314

MAPE is 1.38. This means the model is stable with an accuracy of 1.38%.

Testing the model

acf(tsdataholtfullforecast$residuals,na.action = na.pass,lag.max = 20)

The correlogram shows that the autocorrelations for the in-sample forecast errors do not exceed thesignificance bounds for lags 1-20.

 Box.test(tsdataholtfullforecast$residuals, lag=20, type="Ljung-Box")

    Box-Ljung test

data:  tsdataholtfullforecast$residuals
X-squared = 17.53, df = 20, p-value = 0.6183

Here the Ljung-Box test statistic is 17.53, and the p-value is 0.6183, so there is little evidence of non-zero autocorrelations in the in-sample forecast errors at lags 1-20.

To be sure that the predictive model cannot be improved upon, it is also a good idea to check whether the forecast errors are normally distributed with mean zero and constant variance. To check whether the forecast errors have constant variance, we can make a time plot of the in-sample forecast errors:

We can check whether the forecast errors have constant variance over time, and are normally distributed withmean zero, by making a time plot of the forecast errors

plot(tsdataholtfullforecast$residuals)

The plot shows that the in-sample forecast errors seem to have roughly constant variance over time.

checkresiduals(tsdataholtfullforecast)

pacf(tsdataholtfullforecast$residuals,na.action = na.pass,lag.max = 20)

Arima Model

tsdataar <- ts(data,start=c(1987,1),end=c(1993,12),frequency = 12)

tsdataar

           Jan       Feb       Mar       Apr       May       Jun       Jul
1987   1664.81   2397.53   2840.71   3547.29   3752.96   3714.74   4349.61
1988   2499.81   5198.24   7225.14   4806.03   5900.88   4951.34   6179.12
1989   4717.02   5702.63   9957.58   5304.78   6492.43   6630.80   7349.62
1990   5921.10   5814.58  12421.25   6369.77   7609.12   7224.75   8121.22
1991   4826.64   6470.23   9638.77   8821.17   8722.37  10209.48  11276.55
1992   7615.03   9849.69  14558.40  11587.33   9332.56  13082.09  16732.78
1993  10243.24  11266.88  21826.84  17357.33  15997.79  18601.53  26155.15
           Aug       Sep       Oct       Nov       Dec
1987   3566.34   5021.82   6423.48   7600.60  19756.21
1988   4752.15   5496.43   5835.10  12600.08  28541.72
1989   8176.62   8573.17   9690.50  15151.84  34061.01
1990   7979.25   8093.06   8476.70  17914.66  30114.41
1991  12552.22  11637.39  13606.89  21822.11  45060.69
1992  19888.61  23933.38  25391.35  36024.80  80721.71
1993  28586.52  30505.41  30821.33  46634.38 104660.67

plot(tsdataar)

abline(reg=lm(tsdata~time(tsdata)))

Time series stationarity check - KPSS Test

kpss.test(tsdataar)

p-value smaller than printed p-value

    KPSS Test for Level Stationarity

data:  tsdataar
KPSS Level = 1.3089, Truncation lag parameter = 3, p-value = 0.01

The pvalue from the kpss.test is 0.01, which tells that the series is non-stationary

Making the variance constant by taking log

plot(log(tsdataar))

abline(reg=lm(log(tsdata)~time(tsdata)))

plot(diff(log(tsdataar)))

kpss.test(diff(log(tsdataar)))

p-value greater than printed p-value

    KPSS Test for Level Stationarity

data:  diff(log(tsdataar))
KPSS Level = 0.062948, Truncation lag parameter = 3, p-value = 0.1

The pvalue from the kpss.test is 0.1, which tells that the series is stationary

ACF Plot

acf(diff(log(tsdataar)))

q Value is “2”

PACF Plot

pacf(diff(log(tsdataar)))

p value is also “2”

With the values obtained framing Arima Model Coordinates of ARIMA Model is (2,1,2)

Arimamodel <- arima(log(tsdataar),c(2,1,2),seasonal = list(order=c(2,1,2),period=12))

possible convergence problem: optim gave code = 1

Arimamodel

Call:
arima(x = log(tsdataar), order = c(2, 1, 2), seasonal = list(order = c(2, 1, 
    2), period = 12))

Coefficients:
          ar1      ar2     ma1     ma2     sar1    sar2    sma1     sma2
      -1.2385  -0.5913  0.7873  0.3246  -0.7545  0.1805  0.2113  -0.6426
s.e.   0.3201   0.2813  0.3670  0.3301   0.8286  0.3482  0.8757   0.6065

sigma^2 estimated as 0.02963:  log likelihood = 21.13,  aic = -24.26

Accuracy MAPE ratio for Arima Model

accuracy(Arimamodel)

                       ME      RMSE      MAE         MPE     MAPE      MASE
Training set -0.005248748 0.1582892 0.115803 -0.07967451 1.257853 0.2872351
                    ACF1
Training set -0.01883997

MAPE is 1.25. This means the model is stable with an accuracy of 1.25%

Predicting for next five years

prediction<- predict(Arimamodel,n.ahead = 60)

prediction

$pred
           Jan       Feb       Mar       Apr       May       Jun       Jul
1994  9.604674  9.812547 10.358162 10.073706 10.039623 10.165208 10.421712
1995  9.923647 10.126206 10.650680 10.398505 10.373022 10.502098 10.735892
1996 10.255060 10.472816 10.997617 10.710005 10.684855 10.809559 11.042827
1997 10.562546 10.767821 11.288417 11.033577 11.009554 11.138279 11.367834
1998 10.890228 11.107677 11.631492 11.345537 11.320724 11.445622 11.677888
           Aug       Sep       Oct       Nov       Dec
1994 10.476713 10.565725 10.611994 11.088262 11.853728
1995 10.799250 10.880864 10.950680 11.392429 12.172822
1996 11.095134 11.186420 11.244459 11.724234 12.484964
1997 11.429946 11.512657 11.583785 12.028663 12.806920
1998 11.730605 11.821532 11.880659 12.358729 13.120214

$se
           Jan       Feb       Mar       Apr       May       Jun       Jul
1994 0.1727393 0.1970255 0.2447493 0.2765732 0.3012012 0.3324753 0.3533937
1995 0.5022823 0.5309141 0.5653173 0.5958604 0.6231735 0.6524359 0.6777117
1996 0.8381071 0.8680724 0.9023325 0.9337944 0.9630975 0.9937637 1.0214603
1997 1.1925651 1.2232781 1.2570774 1.2888719 1.3190626 1.3502294 1.3791640
1998 1.5588112 1.5914000 1.6272260 1.6611096 1.6934964 1.7269431 1.7581889
           Aug       Sep       Oct       Nov       Dec
1994 0.3776207 0.3990540 0.4185269 0.4388849 0.4567242
1995 0.7036308 0.7282559 0.7516550 0.7750675 0.7971495
1996 1.0496508 1.0767639 1.1029242 1.1290311 1.1540174
1997 1.4084225 1.4368312 1.4644709 1.4920178 1.5186649
1998 1.7898452 1.8206731 1.8507729 1.8808275 1.9099754

plot(prediction$pred)

Conclusion :

1.Forecast of next Five years has been projected 2.As per both the models performance the ARIMA MODEL FORECAST seems to be better with a MAPE score of 1.25 than Holt Winter’s score of 1.38