Analisis Peramalan Jumlah Penumpang Kereta Api di Jabodetabek dengan Metode Seasonal Autoregressive Intergrated Moving Average (SARIMA)
Novika Sijabat
Maret 23, 2021
Penelitian ini dapat digolongkan kedalam jenis penelitian kuantitatif. Penelitian ini bertujuan untuk mengolah data kuantitatif (angka) yang berupa data jumlah penumpang kereta api periode 2016-2019 di Jabodetabek untuk menghasilkan ramalan jumlah penumpang kereta api 2 tahun mendatang.
Sumber data penelitian ini berasal dari data sekunder. Data sekunder adalah data yang mengacu pada informasi yang dikumpulkan dari sumber yang telah ada. Sumber data sekunder adalah catatan atau dokumentasi perusahaan, publikasi pemerintah, analisis industri oleh media, situs Web, internet dan seterusnya (Uma Sekaran, 2011). a. Studi Kepustakaan, yaitu mencari data buku atau sumber data tertulis lainnya untuk mendapatkan konsep teori dan data yang dapat mendukung analisis dan pembahasan penelitian ini. b. Data yang dipublikasi oleh badan pusat statistik (BPS) yang diperoleh dari PT Kereta Api Indonesia dan PT. KAI Commuter Jabodetabek.
installpackages
library(readr)## Warning: package 'readr' was built under R version 3.6.3library(tsoutliers)## Warning: package 'tsoutliers' was built under R version 3.6.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(forecast)## Warning: package 'forecast' was built under R version 3.6.3library(tseries)## Warning: package 'tseries' was built under R version 3.6.3library(lmtest)## Warning: package 'lmtest' was built under R version 3.6.3## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(FitAR)## Warning: package 'FitAR' was built under R version 3.6.3## Loading required package: lattice## Warning: package 'lattice' was built under R version 3.6.3## Loading required package: leaps## Warning: package 'leaps' was built under R version 3.6.3## Loading required package: ltsa## Loading required package: bestglm## Warning: package 'bestglm' was built under R version 3.6.3##
## Attaching package: 'FitAR'## The following object is masked from 'package:forecast':
##
## BoxCoxinput data
##INPUT DATA
data=read.csv("datapenumpang.csv", header = T, sep = ";")
head(data)## jumlah.penumpang
## 1 22238
## 2 21229
## 3 23206
## 4 23149
## 5 24401
## 6 23821##UBAH DATA MENJADI FORMAT TIME SERIES (+ PLOT DATA)
Yt <- ts(data, start = c(2016,1), freq=12)
Yt## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2016 22238 21229 23206 23149 24401 23821 21574 23923 23570 24533 24104 24841
## 2017 24185 21743 25775 25411 27385 24432 27016 27679 26158 28765 28246 29059
## 2018 28075 25362 29223 28942 28995 24833 29086 28098 27618 29317 28049 29201
## 2019 27768 25305 28366 28062 28369 25816 29714 27651 28293 29278 28563 28860Menguji kestasioneran data dengan melihat plot data (time series plot) dan plot ACF serta PACF.
plot(Yt)
Berdasarkan plot data diatas, plot time series dari penumpang kereta api di Jabodetabek dapat diketahui secara visual bahwa data penumpang kereta api di Jabodetabek dari bulan Januari 2016 hingga Desember 2019 belum stasioner. Namun, kestasioneran data tidak dapat dapat ditentukan hanya dengan menggunakan plot. Oleh karena itu, perlu dilakukan uji kestasioneran data yang meliputi uji kestasioneran dalam varians dan juga uji kestasioneran dalam mean. Uji kestasioneran dalam varians dilakukan dengan melihat nilai lambda menggunakan metode Box-Cox transformation. Sedangkan uji kestasioneran data dalam mean dilakukan dengan menggunakan Augmented Dickey Fuller (ADF) test.
##PLOT ACF PACF
#ACF PACF Jumlah Penumpang
par(mfrow=c(2,1))
acf(Yt, lag=100)
pacf(Yt, lag=100)
##UJI STASIONER TERHADAP VARIANS (BOX.COX)
lambda=BoxCox.lambda(Yt)
lambda## [1] 0.04704793Apabila data yang digunakan belum stasioner dalam varian maka perlu dilakukan transformasi data yaitu dengan menggunakan metode Box-Cox
##TRANSFORMASI DATA
Ytr=((Yt^lambda)-1)/lambda
##UJI STASIONER TERHADAP VARIANS yang sudah ditransformasi(BOX.COX)
lambda=BoxCox.lambda(Ytr)
lambda## [1] 0.9849505par(mfrow=c(2,1))
acf(Ytr, lag=100)
pacf(Ytr, lag=100)
Kemudian apa bila data yang digunakan belum stasioner dalam rata-rata maka perlu dilakukan pembedaan (differencing).
##UJI STASIONER TERHADAP RATA-RATA (ADF)
adf.test(Ytr)##
## Augmented Dickey-Fuller Test
##
## data: Ytr
## Dickey-Fuller = -1.8698, Lag order = 3, p-value = 0.6261
## alternative hypothesis: stationary##DIFFERENCING
Yt1=diff(Ytr, lag=12, differences = 1)
adf.test(Yt1)##
## Augmented Dickey-Fuller Test
##
## data: Yt1
## Dickey-Fuller = -1.9455, Lag order = 3, p-value = 0.5944
## alternative hypothesis: stationaryacf(Yt1, lag=100)
pacf(Yt1, lag=100)
Yt2=diff(Yt1, differences = 1)
adf.test(Yt2)##
## Augmented Dickey-Fuller Test
##
## data: Yt2
## Dickey-Fuller = -2.7772, Lag order = 3, p-value = 0.2713
## alternative hypothesis: stationarypar(mfrow=c(2,1))
acf(Yt2, lag=100)
pacf(Yt2, lag=100)
Setelah memperoleh data yang stasioner dalam varians dan mean, selanjutnya melakukan pendugaan model awal menggunakan plot ACF dan PACF. Plot ACF dan PACF data penumpang kereta api di Jabodetabek setelah differencing dapat di lihat di gambar berikut
Yt3=diff(Yt2, differences = 1)
adf.test(Yt3)## Warning in adf.test(Yt3): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: Yt3
## Dickey-Fuller = -6.2444, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationarypar(mfrow=c(2,1))
acf(Yt3, lag=100)
pacf(Yt3, lag=100)
Lalu tentukan orde pada model dengan melihat plot ACF serta PACF dan menetapkan model dugaan sementara berdasarkan orde model yang telah ditentukan. Setelah itu menguji signifikansi parameter modelnya. Uji signifikansi parameter menggunakan uji t dan p-value dengan kriteria pengambilan keputusan sebagai berikut. Tolak H0 jika │thitung│ > 𝑡𝑎/2 ,df dengan df = T-p dengan T adalah banyaknya data serta p adalah banyaknya parameter ada pada model.
##PENAKSIRAN PARAMETER (4,1,1)(0,0,1)12
fitS1=Arima(Ytr, order=c(0,1,0),seasonal=list(order=c(0,0,1),period=12))
fitS2=Arima(Ytr, order=c(1,1,0),seasonal=list(order=c(0,0,1),period=12))
fitS3=Arima(Ytr, order=c(2,1,0),seasonal=list(order=c(0,0,1),period=12))
fitS4=Arima(Ytr, order=c(3,1,0),seasonal=list(order=c(0,0,1),period=12))
fitS5=Arima(Ytr, order=c(4,1,0),seasonal=list(order=c(0,0,1),period=12))
fitS6=Arima(Ytr, order=c(0,1,1),seasonal=list(order=c(0,0,1),period=12))
fitS7=Arima(Ytr, order=c(1,1,1),seasonal=list(order=c(0,0,1),period=12))
fitS8=Arima(Ytr, order=c(2,1,1),seasonal=list(order=c(0,0,1),period=12))
fitS9=Arima(Ytr, order=c(3,1,1),seasonal=list(order=c(0,0,1),period=12))
fitS10=Arima(Ytr, order=c(4,1,1),seasonal=list(order=c(0,0,1),period=12))
list(fitS1,fitS2,fitS3,fitS4,fitS5,fitS6,fitS7,fitS8,fitS9,fitS10)## [[1]]
## Series: Ytr
## ARIMA(0,1,0)(0,0,1)[12]
##
## Coefficients:
## sma1
## 1.000
## s.e. 0.288
##
## sigma^2 estimated as 0.00564: log likelihood=45.97
## AIC=-87.93 AICc=-87.66 BIC=-84.23
##
## [[2]]
## Series: Ytr
## ARIMA(1,1,0)(0,0,1)[12]
##
## Coefficients:
## ar1 sma1
## -0.5002 0.9995
## s.e. 0.1234 0.3514
##
## sigma^2 estimated as 0.004268: log likelihood=52.91
## AIC=-99.82 AICc=-99.26 BIC=-94.27
##
## [[3]]
## Series: Ytr
## ARIMA(2,1,0)(0,0,1)[12]
##
## Coefficients:
## ar1 ar2 sma1
## -0.6303 -0.2495 1.0000
## s.e. 0.1397 0.1378 0.4461
##
## sigma^2 estimated as 0.004068: log likelihood=54.48
## AIC=-100.96 AICc=-100.01 BIC=-93.56
##
## [[4]]
## Series: Ytr
## ARIMA(3,1,0)(0,0,1)[12]
##
## Coefficients:
## ar1 ar2 ar3 sma1
## -0.6805 -0.3726 -0.1871 0.9999
## s.e. 0.1420 0.1630 0.1376 0.6695
##
## sigma^2 estimated as 0.003999: log likelihood=55.37
## AIC=-100.75 AICc=-99.28 BIC=-91.5
##
## [[5]]
## Series: Ytr
## ARIMA(4,1,0)(0,0,1)[12]
##
## Coefficients:
## ar1 ar2 ar3 ar4 sma1
## -0.6886 -0.3882 -0.2159 -0.0394 0.9999
## s.e. 0.1449 0.1724 0.1728 0.1417 0.9091
##
## sigma^2 estimated as 0.004087: log likelihood=55.41
## AIC=-98.82 AICc=-96.72 BIC=-87.72
##
## [[6]]
## Series: Ytr
## ARIMA(0,1,1)(0,0,1)[12]
##
## Coefficients:
## ma1 sma1
## -0.5887 0.9999
## s.e. 0.0992 0.5862
##
## sigma^2 estimated as 0.00395: log likelihood=54.65
## AIC=-103.3 AICc=-102.74 BIC=-97.75
##
## [[7]]
## Series: Ytr
## ARIMA(1,1,1)(0,0,1)[12]
##
## Coefficients:
## ar1 ma1 sma1
## -0.1840 -0.4863 0.999
## s.e. 0.2046 0.1684 0.639
##
## sigma^2 estimated as 0.003974: log likelihood=55.03
## AIC=-102.06 AICc=-101.1 BIC=-94.65
##
## [[8]]
## Series: Ytr
## ARIMA(2,1,1)(0,0,1)[12]
##
## Coefficients:
## ar1 ar2 ma1 sma1
## -0.2762 -0.0879 -0.4018 0.9997
## s.e. 0.3115 0.2162 0.2849 0.6744
##
## sigma^2 estimated as 0.004049: log likelihood=55.11
## AIC=-100.21 AICc=-98.75 BIC=-90.96
##
## [[9]]
## Series: Ytr
## ARIMA(3,1,1)(0,0,1)[12]
##
## Coefficients:
## ar1 ar2 ar3 ma1 sma1
## -0.5976 -0.3208 -0.1676 -0.0865 1.0000
## s.e. 0.4948 0.3414 0.1855 0.4909 0.7452
##
## sigma^2 estimated as 0.004091: log likelihood=55.39
## AIC=-98.78 AICc=-96.68 BIC=-87.68
##
## [[10]]
## Series: Ytr
## ARIMA(4,1,1)(0,0,1)[12]
##
## Coefficients:
## ar1 ar2 ar3 ar4 ma1 sma1
## -1.3784 -0.8813 -0.5322 -0.2463 0.6868 0.8660
## s.e. 0.2890 0.2889 0.2519 0.1480 0.2601 0.6362
##
## sigma^2 estimated as 0.004623: log likelihood=55.96
## AIC=-97.92 AICc=-95.05 BIC=-84.97##UJI SIGNIFIKANSI PARAMETER
printstatarima=function(Ytr,digits=4,se=T,...){
if (length(Ytr$coef)>0){
cat("\nCoefficients:\n")
coef=round(Ytr$coef,digits=digits)
if(se && nrow(Ytr$var.coef)){
ses=rep(0, length(coef))
ses[Ytr$mask] = round(sqrt(diag(Ytr$var.coef)),digits=digits)
coef=matrix(coef,1,dimnames=list(NULL, names(coef)))
coef=rbind(coef,s.e.=ses)
statt=coef[1,]/ses
pval=2*pt(abs(statt),df=length(Ytr$residuals)-1,lower.tail=F)
coef=rbind(coef,t=round(statt,digits=digits),sign.=round(pval,digits=digits))
coef=t(coef)
}
print.default(coef,print.gap=2)
}
}
printstatarima(fitS1)##
## Coefficients:
## s.e. t sign.
## sma1 1 0.288 3.4722 0.0011printstatarima(fitS2)##
## Coefficients:
## s.e. t sign.
## ar1 -0.5002 0.1234 -4.0535 0.0002
## sma1 0.9995 0.3514 2.8443 0.0066printstatarima(fitS3)##
## Coefficients:
## s.e. t sign.
## ar1 -0.6303 0.1397 -4.5118 0.0000
## ar2 -0.2495 0.1378 -1.8106 0.0766
## sma1 1.0000 0.4461 2.2416 0.0297printstatarima(fitS4)##
## Coefficients:
## s.e. t sign.
## ar1 -0.6805 0.1420 -4.7923 0.0000
## ar2 -0.3726 0.1630 -2.2859 0.0268
## ar3 -0.1871 0.1376 -1.3597 0.1804
## sma1 0.9999 0.6695 1.4935 0.1420printstatarima(fitS5)##
## Coefficients:
## s.e. t sign.
## ar1 -0.6886 0.1449 -4.7522 0.0000
## ar2 -0.3882 0.1724 -2.2517 0.0291
## ar3 -0.2159 0.1728 -1.2494 0.2177
## ar4 -0.0394 0.1417 -0.2781 0.7822
## sma1 0.9999 0.9091 1.0999 0.2770printstatarima(fitS6)##
## Coefficients:
## s.e. t sign.
## ma1 -0.5887 0.0992 -5.9345 0.0000
## sma1 0.9999 0.5862 1.7057 0.0947printstatarima(fitS7)##
## Coefficients:
## s.e. t sign.
## ar1 -0.1840 0.2046 -0.8993 0.3731
## ma1 -0.4863 0.1684 -2.8878 0.0058
## sma1 0.9990 0.6390 1.5634 0.1247printstatarima(fitS8)##
## Coefficients:
## s.e. t sign.
## ar1 -0.2762 0.3115 -0.8867 0.3798
## ar2 -0.0879 0.2162 -0.4066 0.6862
## ma1 -0.4018 0.2849 -1.4103 0.1650
## sma1 0.9997 0.6744 1.4824 0.1449printstatarima(fitS9)##
## Coefficients:
## s.e. t sign.
## ar1 -0.5976 0.4948 -1.2078 0.2332
## ar2 -0.3208 0.3414 -0.9397 0.3522
## ar3 -0.1676 0.1855 -0.9035 0.3709
## ma1 -0.0865 0.4909 -0.1762 0.8609
## sma1 1.0000 0.7452 1.3419 0.1861printstatarima(fitS10)##
## Coefficients:
## s.e. t sign.
## ar1 -1.3784 0.2890 -4.7696 0.0000
## ar2 -0.8813 0.2889 -3.0505 0.0037
## ar3 -0.5322 0.2519 -2.1127 0.0400
## ar4 -0.2463 0.1480 -1.6642 0.1027
## ma1 0.6868 0.2601 2.6405 0.0112
## sma1 0.8660 0.6362 1.3612 0.1799Diketahui bahwa hanya model SARIMA (0,1,0)(0,0,1) 12 dan SARIMA (1,1,0)(0,0,1) 12 yang parameternya signifikan (P-value berada dibawah level toleransi (α = 0,05). Dengan demikian model tersebut memenuhi syarat signifikansi parameter.
Lakukan diagnostic checking
##UJI ASUMSI RESIDUAL (01)
#Pengujian Non Autokorelasi
resY1=residuals(fitS1)
wnY1=Box.test(resY1, type=c("Ljung-Box"))
wnY1##
## Box-Ljung test
##
## data: resY1
## X-squared = 13.847, df = 1, p-value = 0.0001983#Pengujian Normalitas
nY1=length(resY1)
sdY1=sd(resY1)
resnY1=rnorm(nY1,0,sdY1)
ks.test(resY1,resnY1)##
## Two-sample Kolmogorov-Smirnov test
##
## data: resY1 and resnY1
## D = 0.083333, p-value = 0.9969
## alternative hypothesis: two-sided#Pengujian Homoskedastisitas
Box.test(resY1^2, type=c("Ljung-Box"))##
## Box-Ljung test
##
## data: resY1^2
## X-squared = 1.3224, df = 1, p-value = 0.2502##UJI ASUMSI RESIDUAL (02)
#Pengujian Non Autokorelasi
resY2=residuals(fitS2)
wnY2=Box.test(resY2, type=c("Ljung-Box"))
wnY2##
## Box-Ljung test
##
## data: resY2
## X-squared = 1.1623, df = 1, p-value = 0.281#Pengujian Normalitas
nY2=length(resY2)
sdY2=sd(resY2)
resnY2=rnorm(nY2,0,sdY2)
ks.test(resY2,resnY2)##
## Two-sample Kolmogorov-Smirnov test
##
## data: resY2 and resnY2
## D = 0.1875, p-value = 0.3707
## alternative hypothesis: two-sided#Pengujian Homoskedastisitas
Box.test(resY2^2, type=c("Ljung-Box"))##
## Box-Ljung test
##
## data: resY2^2
## X-squared = 0.15689, df = 1, p-value = 0.692dekomposisi<-decompose(Yt,"multiplicative")
dekomposisi## $x
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2016 22238 21229 23206 23149 24401 23821 21574 23923 23570 24533 24104 24841
## 2017 24185 21743 25775 25411 27385 24432 27016 27679 26158 28765 28246 29059
## 2018 28075 25362 29223 28942 28995 24833 29086 28098 27618 29317 28049 29201
## 2019 27768 25305 28366 28062 28369 25816 29714 27651 28293 29278 28563 28860
##
## $seasonal
## Jan Feb Mar Apr May Jun Jul
## 2016 0.9997982 0.8986086 1.0309072 1.0139561 1.0388996 0.9163046 0.9950453
## 2017 0.9997982 0.8986086 1.0309072 1.0139561 1.0388996 0.9163046 0.9950453
## 2018 0.9997982 0.8986086 1.0309072 1.0139561 1.0388996 0.9163046 0.9950453
## 2019 0.9997982 0.8986086 1.0309072 1.0139561 1.0388996 0.9163046 0.9950453
## Aug Sep Oct Nov Dec
## 2016 1.0197282 0.9851088 1.0453747 1.0128798 1.0433890
## 2017 1.0197282 0.9851088 1.0453747 1.0128798 1.0433890
## 2018 1.0197282 0.9851088 1.0453747 1.0128798 1.0433890
## 2019 1.0197282 0.9851088 1.0453747 1.0128798 1.0433890
##
## $trend
## Jan Feb Mar Apr May Jun Jul Aug
## 2016 NA NA NA NA NA NA 23463.54 23566.08
## 2017 24516.42 24899.67 25164.00 25448.17 25797.08 26145.42 26483.25 26796.12
## 2018 27782.33 27886.04 27964.33 28048.17 28062.96 28060.67 28053.79 28038.62
## 2019 27947.42 27954.96 27964.46 27990.96 28010.75 28017.96 NA NA
## Sep Oct Nov Dec
## 2016 23694.54 23895.83 24114.42 24264.21
## 2017 27090.58 27381.37 27595.58 27679.37
## 2018 28000.54 27928.17 27865.42 27880.29
## 2019 NA NA NA NA
##
## $random
## Jan Feb Mar Apr May Jun Jul
## 2016 NA NA NA NA NA NA 0.9240474
## 2017 0.9866809 0.9717518 0.9935722 0.9847956 1.0218063 1.0198201 1.0251960
## 2018 1.0107382 1.0121061 1.0136797 1.0176652 0.9945259 0.9658092 1.0419565
## 2019 0.9937807 1.0073420 0.9839480 0.9887391 0.9748678 1.0055706 NA
## Aug Sep Oct Nov Dec
## 2016 0.9955058 1.0097807 0.9821018 0.9868575 0.9811981
## 2017 1.0129639 0.9801713 1.0049331 1.0105539 1.0061857
## 2018 0.9827301 1.0012479 1.0041651 0.9937885 1.0038161
## 2019 NA NA NA NA NA
##
## $figure
## [1] 0.9997982 0.8986086 1.0309072 1.0139561 1.0388996 0.9163046 0.9950453
## [8] 1.0197282 0.9851088 1.0453747 1.0128798 1.0433890
##
## $type
## [1] "multiplicative"
##
## attr(,"class")
## [1] "decomposed.ts"autoplot(dekomposisi)
##HASIL (CEK MAPE) untuk model2 yang terpenuhi
summary(fitS2)## Series: Ytr
## ARIMA(1,1,0)(0,0,1)[12]
##
## Coefficients:
## ar1 sma1
## -0.5002 0.9995
## s.e. 0.1234 0.3514
##
## sigma^2 estimated as 0.004268: log likelihood=52.91
## AIC=-99.82 AICc=-99.26 BIC=-94.27
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.007517187 0.0632533 0.05146673 0.05567973 0.395157 0.472785
## ACF1
## Training set -0.1508691Pemilihan Model Terbaik Setelah melalui proses estimasi parameter dan diagnostic checking dapat diketahui bahwa model SARIMA (1,1,0)(0,0,1) 12 adalah model terbaik yang bisa dipakai dalam meramalkan jumlah data Penumpang Kereta Api di Jabodetabek. Karena model tersebut dianggap model yang layak karena parameter-parameter yang ada di dalamnya telah signifikan serta residual-residualnya telah mengandung asumsi non autokorelasi dan berdistribusi normal.
PLOT Hasil Peramalan Jumlah Penumpang Kereta Api dengan Metode Seasonal Autoregressive Intergrated Moving Average (SARIMA)
##PLOT PERAMALAN
fc=forecast(fitS2,h = 24)
fc## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2020 13.17169 13.07846 13.26491 13.02911 13.31426
## Feb 2020 13.09176 12.98785 13.19566 12.93285 13.25066
## Mar 2020 13.18516 13.06060 13.30972 12.99466 13.37566
## Apr 2020 13.17097 13.03385 13.30808 12.96127 13.38066
## May 2020 13.21530 13.06440 13.36620 12.98452 13.44608
## Jun 2020 13.17738 13.01493 13.33983 12.92893 13.42583
## Jul 2020 13.30291 13.12917 13.47664 13.03721 13.56861
## Aug 2020 13.22570 13.04160 13.40980 12.94415 13.50725
## Sep 2020 13.24040 13.04640 13.43440 12.94371 13.53710
## Oct 2020 13.27540 13.07199 13.47880 12.96431 13.58648
## Nov 2020 13.28005 13.06771 13.49239 12.95530 13.60479
## Dec 2020 13.26335 13.04230 13.48439 12.92528 13.60141
## Jan 2021 13.27170 13.01184 13.53156 12.87428 13.66912
## Feb 2021 13.26752 12.98943 13.54562 12.84221 13.69283
## Mar 2021 13.26961 12.96758 13.57165 12.80769 13.73153
## Apr 2021 13.26857 12.94765 13.58948 12.77777 13.75936
## May 2021 13.26909 12.92881 13.60936 12.74868 13.78949
## Jun 2021 13.26883 12.91097 13.62669 12.72153 13.81612
## Jul 2021 13.26896 12.89399 13.64393 12.69550 13.84242
## Aug 2021 13.26889 12.87773 13.66005 12.67066 13.86712
## Sep 2021 13.26893 12.86213 13.67572 12.64679 13.89106
## Oct 2021 13.26891 12.84711 13.69071 12.62382 13.91400
## Nov 2021 13.26892 12.83260 13.70524 12.60163 13.93621
## Dec 2021 13.26891 12.81856 13.71927 12.58015 13.95767par(mfrow=c(1,1))
plot(fc)
fc1=fitted(fitS2)
lines(fc1,col='blue')