Ekonometrika

Exercise 3

Kontak	: \(\downarrow\)
Email	ali.19arifin@gmail.com
Instagram	https://www.instagram.com/arifin.alicia/
RPubs	https://rpubs.com/aliciaarifin/
Nama	Alicia Arifin
NIM	20214920001
Prodi	Statistika, 2021

Dataset

The dataset consists of historical sales data collected over the past two years, including information on sales revenue, promotional activities, pricing strategies, weather conditions, and holidays.

library(tidyverse)
library(lubridate)

set.seed(123)
n<-100
start_date <- ymd("2022-01-01")
end_date <- ymd("2022-04-10")  # Update end date to have 100 days
dates <- seq(start_date, end_date, by = "day")

# simulate predictor variables
promotional_spending <- runif(n, min = 1000, max = 5000)
price <- rnorm(n, mean = 50, sd = 10)
weather_conditions <- sample(c("sunny", "cloudy", "rainy"), size = n, replace = TRUE)

# simulate sales revenue
sales_trend <- 0.1 * seq(1, n)
seasonal_pattern <- sin(seq(1, n) * 2 * pi / 365 * 7) * 100
sales_noise <- rnorm(n, mean = 0, sd = 100)
sales_revenue <- 1000 + sales_trend + seasonal_pattern + sales_noise

# create dataframe
simulated_data <- tibble(
  Date = dates,
  Promotional_Spending = promotional_spending,
  Price = price,
  Weather_Conditions = weather_conditions,
  Sales_Revenue = sales_revenue
)

# Display the first few rows of the dataset
head(simulated_data)

Objective

The objective of this practicum is to: + Develop a regression model to understand the relationship between sales revenue and various predictors such as promotional spending, pricing, and external factors. + Build a time series model to capture the temporal patterns and trends in sales revenue, accounting for seasonality and other time-related effects. + Evaluate and compare the performance of both models in forecasting future sales revenue

Answer

Variabel Dependen : + Sales_Revenue
Variabel Independen : + Date + Promotional_Spending + Price + Weather COnditions

Regression Model

# correlation check
num_var = simulated_data[,-c(1,4)]
cor(num_var)

##                      Promotional_Spending       Price Sales_Revenue
## Promotional_Spending           1.00000000 -0.02653768    0.16949977
## Price                         -0.02653768  1.00000000   -0.06417281
## Sales_Revenue                  0.16949977 -0.06417281    1.00000000

ketika dilihat dari korelasi antara sales revenue dengan promotional spending dan price, nilai korelasi tidak mencapai 0.5. Nilai korelasinya sangat kecil. Jika dilihat dari korelasi saja, variabel independen tidak memiliki korelasi yang kuat dengan variabel dependen. Sehingga bisa dikatakan kedua variabel independen ini tidak bisa digunakan dalam membuat regresi.

summary(lm(Sales_Revenue~Promotional_Spending+Price, data=simulated_data))

## 
## Call:
## lm(formula = Sales_Revenue ~ Promotional_Spending + Price, data = simulated_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -346.64  -74.89    6.72   83.14  328.96 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          984.32350   79.38830  12.399   <2e-16 ***
## Promotional_Spending   0.01981    0.01179   1.681   0.0961 .  
## Price                 -0.83264    1.39316  -0.598   0.5515    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 133.7 on 97 degrees of freedom
## Multiple R-squared:  0.03229,    Adjusted R-squared:  0.01234 
## F-statistic: 1.619 on 2 and 97 DF,  p-value: 0.2035

Setelah dicoba menggunakan regresi linear, p-value >alpha 0.05, yang berarti variabel independen tidak ada yang berhubungan siknifikan dengan variabel dependen.
Lanjut untuk variabel kategorik dengan menggunakan Anova.

anova1 <- aov(Sales_Revenue~Date+Weather_Conditions, data=simulated_data)
summary(anova1)

##                    Df  Sum Sq Mean Sq F value Pr(>F)  
## Date                1   70677   70677   3.981 0.0488 *
## Weather_Conditions  2   16377    8189   0.461 0.6319  
## Residuals          96 1704236   17752                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dari hasil tabel anova juga dikatakan bahwa variabel date memiliki pengaruh siknifikan terhadap variabel dependen. Sedangkan weather_conditions tidak memiliki pengaruh siknifikan terhadap variabel dependen.
Date merupakan tanggal, apakah jika date dihilangkan dapat membuat weather_conditions menjadi siknifikan? mari kita coba.

anova2 <- aov(Sales_Revenue~Weather_Conditions, data=simulated_data)
summary(anova2)

##                    Df  Sum Sq Mean Sq F value Pr(>F)
## Weather_Conditions  2   17860    8930   0.488  0.615
## Residuals          97 1773430   18283

Dari hasil Anova tetap sama, weather_conditions tidak memiliki pengaruh siknifikan terhadap variabel dependen.

Jadi, jika dilihat dari variabel independen, hanya Date yang memiliki pengaruh siknifikan terhadap variabel dependen, sisanya tidak memiliki pengaruh yang siknifikan. Maka bisa kita katakan bahwa model time-series sepertinya akan lebih cocok dibandingkan dengan model regresi biasa.

Model Building

model1 <- lm(Sales_Revenue~., data=simulated_data)
summary(model1)

## 
## Call:
## lm(formula = Sales_Revenue ~ ., data = simulated_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -356.97  -75.64    2.00   92.15  272.70 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)  
## (Intercept)              1.851e+04  9.203e+03   2.012   0.0471 *
## Date                    -9.188e-01  4.826e-01  -1.904   0.0600 .
## Promotional_Spending     1.668e-02  1.188e-02   1.404   0.1637  
## Price                   -1.047e+00  1.424e+00  -0.735   0.4640  
## Weather_Conditionsrainy -1.002e+01  3.631e+01  -0.276   0.7831  
## Weather_Conditionssunny -2.500e+01  3.214e+01  -0.778   0.4386  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 132.8 on 94 degrees of freedom
## Multiple R-squared:  0.07413,    Adjusted R-squared:  0.02488 
## F-statistic: 1.505 on 5 and 94 DF,  p-value: 0.1956

Dari hasil model1 didapatkan nilai p-value 0.19, lebih dari alpha 0.05, tidak terdapat variabel dependen yang memiliki pengaruh siknifikan terhadap variabel independen. Model yang dibuat hanya menghasilkan Multiple R-squared 7.413% yang berarti variabel independen hanya dapat menjelaskan variansi sebanyak 7.413%. NIlai R-square tersebut sangat kecil, terdapat kurang lebih 92%an variansi lain yang hilang atau faktor lain yang bisa menjelaskan variabel dependen.

model2 <- lm(Sales_Revenue~ Date, data=simulated_data)
summary(model2)

## 
## Call:
## lm(formula = Sales_Revenue ~ Date, data = simulated_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -384.99  -79.06   11.31   82.38  303.69 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 18540.259   8741.085   2.121   0.0364 *
## Date           -0.921      0.459  -2.006   0.0476 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 132.5 on 98 degrees of freedom
## Multiple R-squared:  0.03946,    Adjusted R-squared:  0.02965 
## F-statistic: 4.026 on 1 and 98 DF,  p-value: 0.04757

Dari hasil model2 didapatkan nilai p-value 0.047, kurang dari alpha 0.05, terdapat variabel dependen yang memiliki pengaruh siknifikan terhadap variabel independen. Model yang dibuat hanya menghasilkan Adjusted R-squared 0.02965 yang berarti variabel independen hanya dapat menjelaskan variansi sebanyak 2.96%. NIlai R-square tersebut sangat kecil, terdapat kurang lebih 97.04%an variansi lain yang hilang atau faktor lain yang bisa menjelaskan variabel dependen.

library(forecast)
accuracy(model1)

##                        ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 5.691975e-15 128.7827 103.2155 -1.793453 10.79428 0.9794119

accuracy(model2)

##                        ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 1.363396e-14 131.1721 104.0721 -1.863197 10.91572 0.9875402

Dari hasil kedua model regresi yang telah dibuat, regrsi linear terbaik memiliki nilai MAPE sebesar 10.79 (model2).

Time Series Model

Time-series Plot

library(plotly)
plot_ly(simulated_data, type = 'scatter', mode = 'lines')%>%
  add_trace(x = ~Date, y = ~(Sales_Revenue), name = 'Sales')%>%
  layout(showlegend = F)

options(warn = -1)

Data terliat memiliki perubahan terhadap trend dan level. Tidak stasioner terhadap rata-rata, tetapi stasioner terhadap varians. Jika data timeseries tidak stasioner terhadap rata-rata, maka diperlukan differencing.

Plot ACF and PACF

y <- simulated_data$Sales_Revenue
print(acf(y,lag = length(y) - 1,pl=TRUE))

## 
## Autocorrelations of series 'y', by lag
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  1.000  0.294  0.242  0.293  0.352  0.183  0.324  0.274  0.187  0.157  0.173 
##     11     12     13     14     15     16     17     18     19     20     21 
##  0.016  0.067  0.115 -0.066 -0.143 -0.141 -0.080 -0.115 -0.097 -0.240 -0.282 
##     22     23     24     25     26     27     28     29     30     31     32 
## -0.282 -0.270 -0.306 -0.166 -0.247 -0.312 -0.327 -0.290 -0.240 -0.140 -0.213 
##     33     34     35     36     37     38     39     40     41     42     43 
## -0.176 -0.174 -0.167 -0.154 -0.079 -0.021 -0.075  0.037  0.027  0.001  0.067 
##     44     45     46     47     48     49     50     51     52     53     54 
##  0.161  0.034  0.109  0.069  0.177  0.127  0.277  0.141  0.102  0.183  0.162 
##     55     56     57     58     59     60     61     62     63     64     65 
##  0.076  0.189  0.201  0.132  0.098  0.139  0.078  0.076  0.074  0.032 -0.024 
##     66     67     68     69     70     71     72     73     74     75     76 
##  0.056  0.089 -0.026  0.006  0.003 -0.016 -0.073 -0.017 -0.080 -0.066 -0.077 
##     77     78     79     80     81     82     83     84     85     86     87 
## -0.089 -0.062 -0.072 -0.075 -0.080 -0.088 -0.084 -0.045 -0.086 -0.015 -0.055 
##     88     89     90     91     92     93     94     95     96     97     98 
## -0.028 -0.059  0.017 -0.021  0.001 -0.043 -0.002 -0.026  0.003 -0.008  0.003 
##     99 
## -0.006

pacf(y, lag = length(y) - 1, pl = TRUE)

Jika plot langsung dibuat plot autokorelasi (ACF dan PACF), dapat dilihat bahwa tidak meluruh secara langsung, maka dari itu kita bisa coba melakukan differencing dan plot ulang grafiknya.

# make acf and pacf plot again using differencing data
print(acf(diff(y),lag = length(y) - 1,pl=TRUE))

## 
## Autocorrelations of series 'diff(y)', by lag
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  1.000 -0.451 -0.076  0.009  0.128 -0.214  0.142  0.030 -0.052 -0.024  0.103 
##     11     12     13     14     15     16     17     18     19     20     21 
## -0.137  0.014  0.148 -0.080 -0.041 -0.037  0.052 -0.030  0.131 -0.081 -0.046 
##     22     23     24     25     26     27     28     29     30     31     32 
##  0.015  0.027 -0.113  0.160 -0.013 -0.038 -0.040  0.003 -0.066  0.145 -0.074 
##     33     34     35     36     37     38     39     40     41     42     43 
##  0.031 -0.022  0.005 -0.036  0.014  0.069 -0.132  0.096  0.027 -0.071 -0.021 
##     44     45     46     47     48     49     50     51     52     53     54 
##  0.156 -0.154  0.067 -0.086  0.090 -0.134  0.218 -0.070 -0.107  0.083  0.034 
##     55     56     57     58     59     60     61     62     63     64     65 
## -0.139  0.074  0.050 -0.032 -0.042  0.084 -0.065  0.031  0.013 -0.005 -0.076 
##     66     67     68     69     70     71     72     73     74     75     76 
##  0.026  0.094 -0.091  0.033  0.003  0.000 -0.057  0.093 -0.053  0.018  0.002 
##     77     78     79     80     81     82     83     84     85     86     87 
## -0.038  0.029 -0.001 -0.009  0.021 -0.002 -0.051  0.069 -0.070  0.067 -0.038 
##     88     89     90     91     92     93     94     95     96     97     98 
##  0.031 -0.064  0.067 -0.034  0.047 -0.041  0.027 -0.030  0.009 -0.009  0.006

pacf(diff(y), lag = length(y) - 1, pl = TRUE)

Dari hasil grafik ACF dan PACF, kita dapat menentukan AR, I dan MA. AR adalah Autoregression, I adalah Integrated, dan MA adalah moving average. AR dan MA ditentukan dari grafik ACF dan PACF. I didapat dari berapa kali data akan didifferencing. Nilai AR dan MA ditentukan dari berapa banyak tonjolan diatas critical value. AR berdasarkan ACF dan MA berdasarkan PACF.
Dari data dapat ditentukan bahwa AR-nya 0-4, I-nya 1 dan MA-nya 0-4. dengan I nya 1 karena data tidak stasioner terhadap rata-rata (didifferencing 1x).

Dari hasil grafik didapatkan data tidak memiliki season atau musim.

Box-Jenkins Method

ARIMA(p,d,q)
dari hasil di atas, didapatkan :
AR maksimal 4, I = 1 dan MA maksimal 4. maka kita bisa menggunakan auto.arima() dalam library forecast untuk menentukan model ARIMA manan yang terbaik.
### ARIMA modeling

library(forecast)
model_arima1 <- auto.arima(y, d = 1, max.p = 4, max.q = 4, approximation = T); model_arima1

## Series: y 
## ARIMA(0,1,1) 
## 
## Coefficients:
##           ma1
##       -0.7841
## s.e.   0.0557
## 
## sigma^2 = 15308:  log likelihood = -617.44
## AIC=1238.87   AICc=1239   BIC=1244.06

accuracy(model_arima1)

##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -4.880079 122.4802 96.44003 -1.969922 10.14477 0.7423986
##                     ACF1
## Training set -0.03524384

Box.test(model_arima1$residuals, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  model_arima1$residuals
## X-squared = 0.12798, df = 1, p-value = 0.7205

dari hasil ARIMA, didapatkan model ARIMA yang terbaik dengan ARIMA(0,1,1) dengan keakuratan prediksi MAPE sebesar 10.14% (data tidak displit). lalu, kita akan mencoba apabila kita langsung menggunakan auto.arima() dengan tidak menentukan p,d,dan q.

model_arima2 <- auto.arima(y,approximation = T); model_arima2

## Series: y 
## ARIMA(1,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1      ma1      mean
##       0.9198  -0.7377  1001.994
## s.e.  0.0535   0.0803    36.507
## 
## sigma^2 = 15050:  log likelihood = -621.52
## AIC=1251.04   AICc=1251.47   BIC=1261.47

accuracy(model_arima2)

##                     ME     RMSE      MAE      MPE     MAPE      MASE       ACF1
## Training set -1.929319 120.8257 94.75599 -1.80019 9.983429 0.7294348 -0.0317043

Box.test(model_arima2$residuals, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  model_arima2$residuals
## X-squared = 0.10356, df = 1, p-value = 0.7476

dari hasil ARIMA model kedua, didapatkan ARIMA(1,0,1) dengan MAPE sebesar 9.98%.

Best Box-Jenkins Model

Model yang terbaik dari ARIMA untuk data simulasi ini adalah ARIMA(1,0,1) dengan MAPE senilai 9.98%.

ARCH GARCH Model

Is there ARCH ?

library(dynlm)
# Step 1: Estimate mean equation r = beta + error
sim.mean <- dynlm(Sales_Revenue ~ 1, data = simulated_data)

# Step 2: Retrieve the residuals from the former model and square them
ehatsq <- ts(resid(sim.mean)^2)

# Step 3: regress squared residuals on one-lagged squared residuals
sim.arch <- dynlm(ehatsq ~ L(ehatsq), data = ehatsq)

summary(sim.arch)

## 
## Time series regression with "ts" data:
## Start = 2, End = 100
## 
## Call:
## dynlm(formula = ehatsq ~ L(ehatsq), data = ehatsq)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -20254 -16290  -9416   7700 113252 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 20283.1065  3031.7460   6.690 1.44e-09 ***
## L(ehatsq)      -0.1254     0.1012  -1.239    0.218    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 24390 on 97 degrees of freedom
## Multiple R-squared:  0.01558,    Adjusted R-squared:  0.005435 
## F-statistic: 1.536 on 1 and 97 DF,  p-value: 0.2183

library(FinTS)
sim.archTest <- ArchTest(simulated_data$Sales_Revenue, demean = T)
sim.archTest

## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  simulated_data$Sales_Revenue
## Chi-squared = 15.093, df = 12, p-value = 0.2364

Pada hasil tes di atas untuk mencek apakah ada efek ARCH atau tidak, hasil mengatakan bahwa tidak terdapat ARCH pada data. Maka dari itu, nilai MAPE yang dihasilkan ARCH dan GARCH pasti tidak akan akurat karena tidak cocok.
untuk mengetes apakah benar atau tidak kita akan coba menggunakan ARCH dan GARCH. ### ARCH Model

library(rugarch)
actual <- diff(y)

arch_spec <- ugarchspec(variance.model= list(model="sGARCH", garchOrder = c(1,0)), mean.model = list(armaOrder = c(0,0), include.mean=F))
arch_fit <- ugarchfit(spec=arch_spec, data=diff(y))
fitted_volitality <- sigma(arch_fit)
mape_arch <- mean(abs(actual - fitted_volitality)/actual);mape_arch

## [1] 7.593992

MAPE ARCH(1) adalah 759.3% yang berarti hasilnya sangat tidak akurat, melebihi 100%. Maka ARCH tidak akan digunakan untuk memodelkan data.

GARCH Model

garch_spec <- ugarchspec(variance.model= list(model="sGARCH", garchOrder = c(1,1)), mean.model = list(armaOrder = c(1,1), include.mean=F))
garch_fit <- ugarchfit(spec=garch_spec, data=diff(y))

fitted_volitality <- sigma(garch_fit)
mape_garch <- mean(abs(actual - fitted_volitality)/actual);mape_garch

## [1] 6.385514

MAPE GARCH(1,1) adalah 638.5% yang berarti hasilnya sangat tidak akurat, melebihi 100%. Maka GARCH tidak akan digunakan untuk memodelkan data.

Kesimpulan

Dari hasil semua model di atas, semua model time-series memiliki nilai MAPE lebih kecil dibandingkan dengan regresi linear. Maka dari itu, model yang terbaik adalah metode Box-Jenkins dengan ARIMA(1,0,1). Persentase error prediksinya sebesar 9.98%.
Model yang lebih baik adalah prediksi dengan menggunakan time-series. Meskipun regresi linear memiliki nilai prediksi yang kecil juga, banyak variabel faktor yang lain yang memengaruhi variabel dependen.