##CABRILLAS, JEALLAH CHRISTINE c. and GERALDE, KIMBERLY M. ##A Time Series Analysis of “RETAIL PRICES OF WELL MILLED RICE FROM 1990 to 2020 in Misamis Oriental”
library(readr) #the use of ‘readr’ is to provide a fast and friendly way to read rectangular data like ‘csv’ WMR <- read_csv(“C:/Users/itlab3/Desktop/PROJECT/WMR.csv”) #importing the csv data file to RStudio (Rscript) View(WMR)
class(WMR) #checking the data type of the imported csv data file WMR #the data type is data.frame
#TRANFORMATION (transform data.frame into a time series data) WMRt <- ts(WMR[,1], start =c(1990, 1), frequency=12) #transforming the data set into time series data with frequency 12 since it is a monthly data WMRt #viewing the transformed data class(WMRt) #double checking if the data is now transformed
#DECOMPOSITION (to analyze the trend, season, cycle and irregularities) plot(WMRt, main = “Well Milled Rice”, ylab = “Prices”) #historical plot of Retail Prices of Well Milled Rice
#decomposition of trend plot(decompose(WMRt)\(trend, main= "Trend", xlab= "Year", ylab= "Retail Prices of Well Milled Rice") #decomposition of season plot(decompose(WMRt)\)season,main= “Season”, xlab= “Year”, ylab= “Retail Prices of Well Milled Rice”) #decomposition of cycle plot(WMRt, main = “Cycle”, xlab= “Year”, ylab=“Retail Prices of Well Milled Rice”) #decomposition of random/irregularity plot(decompose(WMRt)$random, main= “Random”, xlab= “Year”, ylab= “Retail Prices of Well Milled Rice”)
##SEASONALITY TEST (we’ll conduct two tests to determine whether the time series data set is seasonal or not) library(seastests) # “Seasonality Tests” An overall test for seasonality of a given time series library(fpp2) #“Forecasting: Principles and Practice” packages library(forecast) # Forecasting Functions for Time Series and Linear models library(aTSA) #Alternative Time Series Analysis: a package to analyze time series data especially the univariate time series
autoplot(WMRt, main = “Well Milled Rice”, ylab = “Prices”)
summary(wo(WMRt))#WO-test gives out a TRUE if the series is seasonal and FALSE otherwise (first test) isSeasonal(WMRt) #second test #The result of the two tests does not identify seasonality in the historical plot
#Since there is no seasonality, the appropriate model to be used in Exponential Smoothing is the Holt’s Model
##EXPONENTIAL SMOOTHING (using Holt Smoothing)
autoplot(WMRt, ylab= “Retail Prices of Well Milled Rice”) # plot the Retail Prices of Well Milled Rice
#Data Partitioning #the Retail Prices of Well Milled Rice of Mis. Or. dataset has only 372 records WMRt_training_set <- WMRt[1:288] #using the first 288 records as our training set WMRt_training_set <- ts(WMRt_training_set, start=1990, frequency=12) #transforming the data set into a time series data with frequency=12 since it is a monthly data WMRt_testing_set <- WMRt[289:372] #using the last 83 records as our test set WMRt_testing_set <- ts(WMRt_testing_set, start=2014, frequency=12) #transforming the data set into a time series data with frequency=12 since it is a monthly data WMRt_training_set autoplot(WMRt_training_set, main= “Training Set”, ylab= “Retail Prices (WMR)”) WMRt_testing_set autoplot(WMRt_testing_set, main= “Testing Set”, ylab= “Retail Prices (WMR)”)
#Model Fitting WMRt_naive_model <- naive(WMRt_training_set,h=length(WMRt_testing_set)) #Naive autoplot(WMRt_naive_model)+autolayer(fitted(WMRt_naive_model), series=“Fitted”)+ylab(“Retail Price (WMR)”) accuracy(WMRt_naive_model) #checking the accuracy of naive model
WMRt_ses_model <- ses(WMRt_training_set, h=length(WMRt_testing_set)) #Simple Exponential Smoothing (SES) autoplot(WMRt_ses_model)+autolayer(fitted(WMRt_ses_model),series=“Fitted”)+ylab(“Retail Price (WMR)”) accuracy(WMRt_ses_model) #checking the accuracy of SES model
WMRt_holt_model <- holt(WMRt_training_set,h=length(WMRt_testing_set))#Holt Smoothing autoplot(WMRt_holt_model)+autolayer(fitted(WMRt_holt_model),series=“Fitted”)+ylab(“Retail Price (WMR)”) accuracy(WMRt_holt_model)
WMRt_holtD_model <- holt(WMRt_training_set,h=length(WMRt_testing_set),damped=TRUE) #Damped Holt Smoothing autoplot(WMRt_holtD_model)+autolayer(fitted(WMRt_holtD_model),series=“Fitted”)+ylab(“Retail Price (WMR)”) accuracy(WMRt_holtD_model)
#Model Selection accuracy(WMRt_naive_model,x=WMRt_testing_set) #accuracy of naive model accuracy(WMRt_ses_model,x=WMRt_testing_set)#accuracy of SES model accuracy(WMRt_holt_model,x=WMRt_testing_set)#accuracy of Holt’s model accuracy(WMRt_holtD_model,x=WMRt)#accuracy of Damped Holt’s model
#Diagnostic checkresiduals(WMRt_holt_model)#checking the residuals of Holt Smoothing
#Final Model for Retail Price of Well Milled Rice from 1990 to 2020 in Misamis Oriental WMRt_holt_model_final <- holt(WMRt, h=7) WMRt_holt_model_final$model WMRt_holt_model_final autoplot(WMRt_holt_model_final)
##ARIMA MODELING (Non-Stationary ARMA)
library(fpp2) library(aTSA) library(forecast)
##Modeling the Retail Prices of Well Milled Rice Using Box-Jenkins Method autoplot(WMRt, main= (“Retail Rrices of WMR”)) #Stationary Test: augmented dickey fuller test #\(H_o\): Non-stationary #\(H_a\): Stationary adf.test(WMRt) #Performs the Augmented Dickey-Fuller test for the null hypothesis of a unit root of a univarate time series x
#DIFFERENCING WMRt_diff <- diff(WMRt,differences = 1) autoplot(WMRt_diff,main=“First Difference of WMR Data”,xlab="") adf.test(WMRt_diff)
#ARMA Order acf(WMRt_diff) #autocorrelation of WMRt_diff pacf(WMRt_diff) #patrial autocorrelation of WMRt_diff
#ARIMA(p,d,q) wmr_model1 <- arima(WMRt, order= c(0,1,0)) wmr_model2 <- arima(WMRt, order=c(1,1,0)) wmr_model3 <- arima(WMRt, order=c(1,1,1))
#Model Selection cat(“AIC of Random Walk”) AIC(wmr_model1) cat(“AIC of ARIMA(1,1,0)”) AIC(wmr_model2) cat(“AIC of ARIMA (1,1,1)”) AIC(wmr_model3)
#Diagnostic checkresiduals(model2) #first choice (based from AIC) is ARIMA(1,1,0)
#Forecasting forecast(model2, lead=5)
#auto.arima()function is used to see the fittest model for the data auto.arima(WMRt) WMRt_auto <-auto.arima(WMRt) WMR_auto_arima <- Arima(WMRt_training_set,order=c(0,1,1), include.drift = T) WMR_forecast <- forecast::forecast(WMR_auto_arima, h=length(WMRt_testing_set))
##to identify which model has lesser error with respect to the training set, however, smaller errors with respect to the training set does not guarantee a fit model AIC(WMRt_auto) accuracy(WMR_forecast,x=WMRt_testing_set)
#Diagnostic of ARIMA(0,1,1) with drift checkresiduals(WMRt_auto)
#FORECASTING for ARIMA(0,1,1) with drift forecast::forecast(WMRt_auto,h=7) autoplot(forecast::forecast(WMRt_auto,h=7))