# 1.LOADING DATA
inflation<-read.csv("~/som.csv")
# 2. SELECTING VARIABLE OF INTEREST
inflation<-inflation$CPI
# 3. DESCRIPTIVE STATISTICS
summary(inflation)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 100.0 124.3 158.9 154.8 184.9 218.2
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 100.0 124.3 158.9 154.8 184.9 218.2
library(AdequacyModel)
descriptive(inflation)
## $mean
## [1] 154.7518
##
## $median
## [1] 158.9
##
## $mode
## [1] 110
##
## $variance
## [1] 1189.735
##
## $Skewness
## [1] 0.07614
##
## $Kurtosis
## [1] -1.32706
##
## $minimum
## [1] 100
##
## $maximum
## [1] 218.23
##
## $n
## [1] 123
## $mean
## [1] 154.7518
##
## $median
## [1] 158.9
##
## $mode
## [1] 110
##
## $variance
## [1] 1189.735
##
## $Skewness
## [1] 0.07614
##
## $Kurtosis
## [1] -1.32706
##
## $minimum
## [1] 100
##
## $maximum
## [1] 218.23
##
## $n
## [1] 123
# 4.CHECKING MISSING VALUE
any(is.na(inflation))
## [1] FALSE
## [1] FALSE
# 5. CHECKING TIME SERIES
class(inflation)
## [1] "numeric"
## [1] "numeric"
is.ts(inflation)
## [1] FALSE
## [1] FALSE
#6.converting to time series
inflationts<-ts(inflation,frequency = 12,start=c(2014,12))
is.ts(inflationts)
## [1] TRUE
## [1] TRUE
#7 ploting
windows()
plot.ts(inflationts,col="pink",
xlab="month",
ylab="inflation",
main="inflation data")
#8. DECOMPOSTION
windows()
add<-plot(decompose(inflationts,type = "additive"))
windows()
mult<-plot(decompose(inflationts,type = "multiplicative"))
#9. CHECKING STATIONARITY
library(tseries)
## Warning: package 'tseries' was built under R version 4.2.3
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
## Warning: package 'tseries' was built under R version 4.2.3
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
adf.test(inflationts)
##
## Augmented Dickey-Fuller Test
##
## data: inflationts
## Dickey-Fuller = -2.9364, Lag order = 4, p-value = 0.1875
## alternative hypothesis: stationary
##
## Augmented Dickey-Fuller Test
##
## data: inflationts
## Dickey-Fuller = -2.9364, Lag order = 4, p-value = 0.1875
## alternative hypothesis: stationary
pp.test(inflationts)
##
## Phillips-Perron Unit Root Test
##
## data: inflationts
## Dickey-Fuller Z(alpha) = -8.39, Truncation lag parameter = 4, p-value =
## 0.6311
## alternative hypothesis: stationary
##
## Phillips-Perron Unit Root Test
##
## data: inflationts
## Dickey-Fuller Z(alpha) = -8.39, Truncation lag parameter = 4, p-value =
## 0.6311
## alternative hypothesis: stationary
# 10.DIFFERENCING
inflationdiff<-diff(inflationts,differences = 1)
#11.CHECKING DIFFERENCING DATA
adf.test(inflationdiff)
##
## Augmented Dickey-Fuller Test
##
## data: inflationdiff
## Dickey-Fuller = -4.0199, Lag order = 4, p-value = 0.01064
## alternative hypothesis: stationary
##
## Augmented Dickey-Fuller Test
##
## data: inflationdiff
## Dickey-Fuller = -4.0199, Lag order = 4, p-value = 0.01064
## alternative hypothesis: stationary
pp.test(inflationdiff)
## Warning in pp.test(inflationdiff): p-value smaller than printed p-value
##
## Phillips-Perron Unit Root Test
##
## data: inflationdiff
## Dickey-Fuller Z(alpha) = -92.094, Truncation lag parameter = 4, p-value
## = 0.01
## alternative hypothesis: stationary
## Warning in pp.test(inflationdiff): p-value smaller than printed p-value
##
## Phillips-Perron Unit Root Test
##
## data: inflationdiff
## Dickey-Fuller Z(alpha) = -92.094, Truncation lag parameter = 4, p-value
## = 0.01
## alternative hypothesis: stationary
#12. GOVERNING NUMBER OF DIFFERENCING
library(forecast)
## Warning: package 'forecast' was built under R version 4.2.3
## Warning: package 'forecast' was built under R version 4.2.3
ndiffs(inflationts)
## [1] 1
## [1] 1
pacf(inflationts)
#13. governing the number of p and q parameters
acf(inflationdiff)
pacf(inflationdiff)
pacf
## function (x, lag.max, plot, na.action, ...)
## UseMethod("pacf")
## <bytecode: 0x000002231172b768>
## <environment: namespace:stats>
## function (x, lag.max, plot, na.action, ...)
## UseMethod("pacf")
## <bytecode: 0x000002b1b425d108>
## <environment: namespace:stats>
# 14. Auto_arima
arimamodel<-auto.arima(inflationts)
summary(arimamodel)
## Series: inflationts
## ARIMA(1,1,0) with drift
##
## Coefficients:
## ar1 drift
## 0.3938 0.9973
## s.e. 0.0902 0.1775
##
## sigma^2 = 1.447: log likelihood = -194.71
## AIC=395.42 AICc=395.62 BIC=403.83
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.01953941 1.188035 0.8304844 -0.0272713 0.5959843 0.07442398
## ACF1
## Training set -0.06762716
## Series: inflationts
## ARIMA(1,1,0) with drift
##
## Coefficients:
## ar1 drift
## 0.3938 0.9973
## s.e. 0.0902 0.1775
##
## sigma^2 = 1.447: log likelihood = -194.71
## AIC=395.42 AICc=395.62 BIC=403.83
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.01953941 1.188035 0.8304844 -0.0272713 0.5959843 0.07442398
## ACF1
## Training set -0.06762716
#15.COMPARING AUTO ARIMA MODELS
arima1<-arima(inflationts,order = c(1,1,2))
## Warning in arima(inflationts, order = c(1, 1, 2)): possible convergence problem:
## optim gave code = 1
## Warning in arima(inflationts, order = c(1, 1, 2)): possible convergence problem:
## optim gave code = 1
arima2<-arima(inflationts,order = c(0,1,2))
arima3<-arima(inflationts,order = c(1,1,0))
arima4<-arima(inflationts,order = c(0,1,0))
aic<-cbind(arima1$aic,arima2$aic,arima3$aic,arima4$aic)
aic
## [,1] [,2] [,3] [,4]
## [1,] 403.4526 426.2482 411.9608 464.1167
## [,1] [,2] [,3] [,4]
## [1,] 403.4526 426.2482 411.9608 464.1167
bic<-cbind(BIC(arima1),BIC(arima2),BIC(arima3),BIC(arima4))
bic
## [,1] [,2] [,3] [,4]
## [1,] 414.6686 434.6602 417.5688 466.9207
## [,1] [,2] [,3] [,4]
## [1,] 414.6686 434.6602 417.5688 466.9207
# 16. forecasting
arima1f<-forecast(arima(inflationts,order = c(1,1,2)))
## Warning in arima(inflationts, order = c(1, 1, 2)): possible convergence problem:
## optim gave code = 1
## Warning in arima(inflationts, order = c(1, 1, 2)): possible convergence problem:
## optim gave code = 1
arima2f<-forecast(arima(inflationts,order = c(0,1,2)))
arima3f<-forecast(arima(inflationts,order = c(1,1,0)))
arima4f<-forecast(arima(inflationts,order = c(0,1,0)))
# forecasting metrics
library(TSstudio)
## Warning: package 'TSstudio' was built under R version 4.2.3
## Warning: package 'TSstudio' was built under R version 4.2.3
splitdata<-ts_split(inflationts,sample.out = 12)
traindata<-splitdata$train
testdata<-splitdata$test
arima1f<-forecast(arima(traindata,order = c(1,1,2)))
## Warning in arima(traindata, order = c(1, 1, 2)): possible convergence problem:
## optim gave code = 1
## Warning in arima(traindata, order = c(1, 1, 2)): possible convergence problem:
## optim gave code = 1
arima2f<-forecast(arima(traindata,order = c(1,1,0)))
arima3f<-forecast(arima(traindata,order = c(0,1,2)))
arima4f<-forecast(arima(traindata,order = c(0,1,0)))
library(Metrics)
## Warning: package 'Metrics' was built under R version 4.2.3
##
## Attaching package: 'Metrics'
## The following object is masked from 'package:forecast':
##
## accuracy
## Warning: package 'Metrics' was built under R version 4.2.3
##
## Attaching package: 'Metrics'
## The following object is masked from 'package:forecast':
##
## accuracy
s1f<-forecast::accuracy(arima1f,testdata)
s2f<-forecast::accuracy(arima2f,testdata)
s3f<-forecast::accuracy(arima3f,testdata)
s4f<-forecast::accuracy(arima4f,testdata)
metrics<-cbind(s1f,s2f,s3f,s4f)
windows()
plot(arima1f)
windows()
plot(arima2f)
windows()
plot(arima3f)
windows()
plot(arima4f)
#comparing auto-arima to other time series models
library(forecast)
arimamodel<-auto.arima(traindata)
arimamodel
## Series: traindata
## ARIMA(2,1,1) with drift
##
## Coefficients:
## ar1 ar2 ma1 drift
## -0.5571 0.2410 0.9803 0.9080
## s.e. 0.1645 0.1546 0.1181 0.1711
##
## sigma^2 = 1.481: log likelihood = -176.19
## AIC=362.39 AICc=362.96 BIC=375.89
## Series: traindata
## ARIMA(2,1,1) with drift
##
## Coefficients:
## ar1 ar2 ma1 drift
## -0.5571 0.2410 0.9803 0.9080
## s.e. 0.1645 0.1546 0.1181 0.1711
##
## sigma^2 = 1.481: log likelihood = -176.19
## AIC=362.39 AICc=362.96 BIC=375.89
etsmodel<-ets(traindata)
etsmodel
## ETS(A,Ad,N)
##
## Call:
## ets(y = traindata)
##
## Smoothing parameters:
## alpha = 0.9999
## beta = 0.181
## phi = 0.9732
##
## Initial states:
## l = 100.3368
## b = 0.1034
##
## sigma: 1.3383
##
## AIC AICc BIC
## 594.3276 595.1353 610.5847
## ETS(A,Ad,N)
##
## Call:
## ets(y = traindata)
##
## Smoothing parameters:
## alpha = 0.9999
## beta = 0.181
## phi = 0.9732
##
## Initial states:
## l = 100.3368
## b = 0.1034
##
## sigma: 1.3383
##
## AIC AICc BIC
## 594.3276 595.1353 610.5847
batsmodel<-bats(traindata)
batsmodel
## BATS(1, {0,0}, 1, -)
##
## Call: bats(y = traindata)
##
## Parameters
## Lambda: 1
## Alpha: 0.9996527
## Beta: 0.06537387
## Damping Parameter: 1
##
## Seed States:
## [,1]
## [1,] 1.071639e+02
## [2,] 3.575319e-03
## attr(,"lambda")
## [1] 0.9999999
##
## Sigma: 1.563269
## AIC: 631.9429
## BATS(1, {0,0}, 1, -)
##
## Call: bats(y = traindata)
##
## Parameters
## Lambda: 1
## Alpha: 0.9996527
## Beta: 0.06537387
## Damping Parameter: 1
##
## Seed States:
## [,1]
## [1,] 1.071639e+02
## [2,] 3.575319e-03
## attr(,"lambda")
## [1] 0.9999999
##
## Sigma: 1.563269
## AIC: 631.9429
tbatsmodel<-tbats(traindata)
tbatsmodel
## BATS(1, {0,0}, 1, -)
##
## Call: tbats(y = traindata)
##
## Parameters
## Alpha: 0.996089
## Beta: 0.03795683
## Damping Parameter: 1
##
## Seed States:
## [,1]
## [1,] 108.7886675
## [2,] 0.3960777
##
## Sigma: 1.592509
## AIC: 634.0568
## BATS(1, {0,0}, 1, -)
##
## Call: tbats(y = traindata)
##
## Parameters
## Alpha: 0.996089
## Beta: 0.03795683
## Damping Parameter: 1
##
## Seed States:
## [,1]
## [1,] 108.7886675
## [2,] 0.3960777
##
## Sigma: 1.592509
## AIC: 634.0568
# FORCASTING FOR SARIMA,ETS,BATS,TBATS
arimaforecast<-forecast(arimamodel,h=12)
etsforecast<-forecast(etsmodel,h=12)
batsforecast<-forecast(batsmodel,h=12)
tbatsforecast<-forecast(tbatsmodel,h=12)
windows()
plot(arimaforecast)
windows()
plot(etsforecast)
windows()
plot(batsforecast)
windows()
plot(tbatsforecast)
