-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathtime-series.Rmd
1881 lines (1317 loc) · 81.5 KB
/
time-series.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
---
title: "Supply Chain Analytics Individual Assignment"
author:
- Candidate Number - 01437771
- Teacher - Dr. Jiahua Wu
- Avi Mago
subtitle: BS1808 Logistics and Supply Chain Analytics
fontsize: 11pt
output:
html_document:
theme: "paper"
number_sections: yes
toc: yes
toc_depth: 4
csl: harvard-imperial-college-london.csl
bibliography: sca_final.bib
---
```{r Loading Libraries, include=FALSE}
library(forecast)
library(tseries)
library(ggplot2)
library(knitr)
library(stargazer)
library(ggthemes)
library(xtable)
library(texreg)
library(kableExtra)
library(gridExtra)
```
```{r Loading dataset, include=FALSE}
# loading transaction data
transaction_rep <- read.csv("results/rep_transactions.csv")
transaction_sub <- read.csv("results/sub_transactions.csv")
# Loading data for stores
data_12631 <- read.csv("results/12631.csv")
data_20974 <- read.csv("results/20974.csv")
data_46673 <- read.csv("results/46673.csv")
data_4904 <- read.csv("results/4904.csv")
```
```{r Analysing existence of outlier, include=FALSE}
# Histogram to see if there are any outliers
repice_hist <- ggplot(data = transaction_rep, aes(x = total_lettuce_trans_rep)) + geom_histogram() +labs(x = "Total lettuce Ussage in Ounces", y = "Number of Such Transaction", title = "Recipe") + theme_few()
sub_hist <- ggplot(data = transaction_sub, aes(x = total_lettuce_trans_sub)) + geom_histogram() + labs(x = "Total lettuce Ussage in Ounces", y = "Number of Such Transaction", title = "Sub-Recipe") + theme_few()
# Outlier Cutoff
# Recipe
rep_cutoff <- mean(transaction_rep$total_lettuce_trans_rep) + (3*sd(transaction_rep$total_lettuce_trans_rep))
# Sub Recipe
sub_cutoff <- mean(transaction_sub$total_lettuce_trans_sub) + (3*sd(transaction_sub$total_lettuce_trans_sub))
# Sub recipe Outlier Pattern check
sub_out <- ggplot(data = transaction_sub, aes(y = total_lettuce_trans_sub, x = date)) + geom_point(alpha = 0.4) + facet_grid(~storenumber) + labs(x = "Date tick marks", y = "Lettuce per Transaction", title = "Checking for patterns in outliers Sub-Recipe")+ geom_hline(aes(yintercept=sub_cutoff), colour="green", linetype="dashed") + theme_few() +theme(axis.text.x = element_blank())
# Recipe Outlier Pattern Check
repice_out <- ggplot(data = transaction_rep, aes(y = total_lettuce_trans_rep, x = date)) + geom_point(alpha = 0.4) + facet_grid(~storenumber) + labs(x = "Date tick marks", y = "Lettuce per Transaction", title = "Checking for patterns in outliers Recipe") + theme(axis.text.x = element_blank())+ geom_hline(aes(yintercept=rep_cutoff), colour="green", linetype="dashed") + theme_few() +theme(axis.text.x = element_blank())
# Plot
# grid.arrange(sub_out, repice_out, nrow=2)
# Green line show the outlier cutoff by using mean + 3 sd formula
# Blue line is choosen by visually accessing the plots and removing the
# Irregular and extermely high value points.
```
```{r Store 12631, include=FALSE}
# Creating time series object
ts_12631 <- ts((data_12631$total_lettuce), frequency = 7, start = c(10, 1))
# Splitting into train and test
ts_12631_train <- subset(ts_12631, end = 93)
ts_12631_test <- subset(ts_12631, start = 94)
```
```{r 12631 Plotting, include=FALSE}
# Plot
# autoplot(ts_12631_train, main = "Store 12631 Time Series Plot", xlab = "Weeks", ylab = "Lettuce Ussage") + theme_economist()
# Seasonality plot
# ggseasonplot(ts_12631_train, main = "Store 12631 Time Series Plot (Seasonality check)", xlab = "Weeks") + theme_economist()
# Sub series plot
# ggsubseriesplot(ts_12631_train , main = "Store 12631 Sub Series Plot", xlab = "Weeks", ylab = "Lettuce Demand") + theme_economist()
```
```{r 12631 ARIMA Stationary Test, include=FALSE}
# stationary test
# Dickey Fuller Test
adf_12631 <- adf.test(log(ts_12631_train))
# Rejected the null that the sereis contain a unit root in favour of stationarity
# Phillips - Perron Test
pp_12631 <- pp.test(log(ts_12631_train))
# Rejected the null that the sereis contain a unit root in favour of stationarity
# KPSS test
kpss_12631 <- kpss.test(log(ts_12631_train))
# Rejected the null that the proces is STATIONARY
# Checking how many differences are needed for the process to be stationary
#since KPSS rejected the null
ndiffs(log(ts_12631_train))
# 1 differences in needed
# Testing again after taking difference
ts_12631_train.diff1 <- diff(log(ts_12631_train), differences = 1)
# Dickey Fuller Test
adf_12631.diff1 <- adf.test(ts_12631_train.diff1)
# Rejected the null that the sereis contain a unit root in favour of stationarity
# Phillips - Perron Test
pp_12631.diff1 <- pp.test(ts_12631_train.diff1)
# Rejected the null that the sereis contain a unit root in favour of stationarity
# KPSS test
kpss_12631.diff1 <- kpss.test(ts_12631_train.diff1)
# Fail to reject the null that the proces is STATIONARY
# Creating a table to summarize first test results
test_results_12631 <- data.frame(Test=c("Dickey Fuller", "Phillips - Perron", "KPSS"),
"P Value" = c(adf_12631$p.value,pp_12631$p.value,kpss_12631$p.value),
"Null Hypothesis" = c("Unit Root", "Unit Root", "Stationary"),
Conclusion = c("Stationary", "Stationary", "Not Stationary"))
# Creating a table to summarize first test results after taking 1st difference
test_results_12631.diff1 <- data.frame(Test=c("Dickey Fuller", "Phillips - Perron", "KPSS"),
"P Value" = c(adf_12631.diff1$p.value,pp_12631.diff1$p.value,
kpss_12631.diff1$p.value),
"Null Hypothesis" = c("Unit Root", "Unit Root", "Stationary"),
Conclusion = c("Stationary", "Stationary", "Stationary"))
# Final test results
final_test_results_12631 <- rbind(test_results_12631, test_results_12631.diff1)
row.names(final_test_results_12631) <- c(1, 2, 3, 1.1, 2.1, 3.1)
# Creating a table to present these results after taking 1st difference
# kable(final_test_results_12631, caption = "Stationary Test Results for Store 12631", format = "html", booktabs = T, align = 'l') %>% kable_styling() %>% group_rows("Before 1st diff", 1, 3, latex_gap_space = "2em") %>% group_rows("After 1st diff", 4, 6, latex_gap_space = "2em")
# Plotting the time series to see if stationary
# autoplot(ts_12631_train.diff1, main = "Store 12631 Time Series Plot (Stationary)", xlab = "Weeks", ylab = "Lettuce Ussage") + theme_economist()
# Looks stationary
```
```{r 12631 ARIMA Plotting time series data, include=FALSE}
# PACF AND ACF
# ggtsdisplay(ts_12631_train.diff1, lag.max = 49, main = "Store 12631 Time Series (Stationary), ACF and PACF Plots", xlab = "Weeks", theme = theme_economist())
# Looking at the ACF and PACF model?
# p = 0, q = 1 (ACF cuts of at 1 and PACF tails off)
# Seasonal - ACF cuts of at 4 and PACF no tailing off Q = 4, P = 0
# Let's see if there is seasonality
nsdiffs(ts_12631_train)
# There is none.
```
```{r 12631 ARIMA Chosing best model, include=FALSE, results='asis'}
auto.arima(ts_12631_train, trace = TRUE, ic = 'bic', lambda = 0)
# Best model is ARIMA(0,1,1)(0,0,1) with drift (BIC - -60.24)
# Using PACF AND ACF Model we get ARIMA(0,1,1)(0,0,4) and ARIMA (0,1,1)(0,0,2)
# Using ARIMA (0,1,1)(0,0,2) we get a BIC -59.25 which is clearly lower than -60.24
# and therefore this is our best model. Let's use both models on the
# test set to check performance
ts_12631_train.arima.m1 <- Arima(ts_12631_train, order = c(0, 1, 1),
seasonal = list(order = c(0, 0, 2), period = 7), lambda = 0)
ts_12631_train.arima.m2 <- Arima(ts_12631_train, order = c(0, 1, 1),
seasonal = list(order = c(0, 0, 1), period = 7), lambda = 0)
# texreg(list(ts_12631_train.arima.m1, ts_12631_train.arima.m2), caption = "ARIMA Self selected V/s automatically selected model for Store 12631",custom.model.names = c("Self", "Auto Arima"))
```
```{r 12631 ARIMA Checking residual, include=FALSE}
# Plot
# Model 1
# checkresiduals(ts_12631_train.arima.m1, xlab = "weeks", theme = theme_economist(), test=FALSE)
# Model 2
# checkresiduals(ts_12631_train.arima.m2, xlab = "Weeks", theme = theme_economist(), test=FALSE)
# For both model the residual looks fine and the p-value for Ljung-Box test is
# way higher than standard 0.05 (m1 - .2764, m2 - .4851)
# Ljung-Box test
ljung_box_test_12631 <- data.frame(Model = c("ARIMA (0,1,1)(0,0,2)", "ARIMA(0,1,1)(0,0,1)"),
"p-value" = c(0.6218, 0.1174))
# Print output
# kable(ljung_box_test_12631, caption = "Ljung-Box test result for the two candidate models (Store 12631)", format = "latex", booktabs = T, align = 'l')
```
```{r 12631 ARIMA accuracy test, include=FALSE}
# Accuracy test for best model
acc_12631_arima.m1 <- accuracy(forecast(ts_12631_train.arima.m1, h = 10), ts_12631_test)
# Accuracy test for best aacording to auto arima test
acc_12631_arima.m2 <- accuracy(forecast(ts_12631_train.arima.m2, h = 10), ts_12631_test)
# Creating a dataframe and renaming the rows to avoid duplicates
acc_12631_arima <- as.data.frame(rbind(acc_12631_arima.m1, acc_12631_arima.m2))
row.names(acc_12631_arima) <- c("Training M1", "Test M1", "Training M2", "Test M2")
# Creating a table for combined accuracy results
# kable(acc_12631_arima, caption = "Out of Sample Performance for Self selected Model V/s Auto Arima Model", format = "latex", booktabs = T) %>% kable_styling() %>% group_rows("Self selected", 1, 2, latex_gap_space = "2em") %>% group_rows("Auto Arima", 3, 4, latex_gap_space = "2em")
# model 1 outperforms the model 2. Therrefore we choose model 1.
```
```{r 12631 ARIMA Training on all data and predicting, include=FALSE, results='asis'}
# Training on all data
ts_12631.arima.m <- Arima(ts_12631, order = c(0, 1, 1),
seasonal = list(order = c(0, 0, 2), period = 7), lambda = 0)
# Forecast for next two weeks
ts_12631.arima.f <- forecast(ts_12631.arima.m, h = 14)
# Final trained model output
# texreg(ts_12631.arima.m, caption = "Final Model for Store 12631 (Trained on all data)")
# Predictions output formatting
ts_12631.arima.f.df <- as.data.frame(ts_12631.arima.f)
days <- data.frame(Day = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14))
final_forecast_12631 <- cbind(days, ts_12631.arima.f.df$`Point Forecast`)
# Final Forecast
# kable(final_forecast_12631, caption = "2 Week Prediction for Store 12631", col.names = c("Day", "Forecast for Store 12631"), align = 'c')
```
```{r 12631 Holt Winters Model, include=FALSE, results='asis'}
# Range bar plots
# ts_12631_train %>% stl(s.window = "period") %>% autoplot + theme_economist() + labs(x = "Weeks")
# Seasonality and trend has a magnitude which is less than the stochastic component.
# Seasonality ranges from +30 to -30. Indicating some days have higher sales than others
# Trend shows a continuos rise till week 18. Followed by a drastic dip in sales in week
# 22 to 24.
# Looking at the range bar plots we can see additive seasonality and linear(additive)
# trend
ts_12631.ets1 <- ets(ts_12631_train, model = "AAA")
# Automatic approach
ts_12631.ets2 <- ets(ts_12631_train, model = "ZZZ")
# According to BIC model 2 outperforms model 1
# Output
# texreg(list(ts_12631.ets1, ts_12631.ets2), caption = "Holt Winters Self selected(A,A,A) V/s automatically selected model(M,N,M) for Store 12631", model.names = c("Self", "Automatic"), single.row = T)
```
```{r 12631 Holt Winters Model Out of sample accuracy, include=FALSE}
# Forecast
# Model 1
ts_12631.ets1.f <- forecast(ts_12631.ets1, h = 10)
# Model 2
ts_12631.ets2.f <- forecast(ts_12631.ets2, h = 10)
# Out of sample accuracy
# Model 1
acc_12631_hw.m1 <-accuracy(ts_12631.ets1.f, ts_12631_test)
# Model 2
acc_12631_hw.m2 <-accuracy(ts_12631.ets2.f, ts_12631_test)
# Model 1 performs better according to out of sample accuracy with most of the accuracy
# measures
# favouring it. Thus we go ahead with Model 1.
# Creating a dataframe and renaming the rows to avoid duplicates
acc_12631_hw <- as.data.frame(rbind(acc_12631_hw.m1, acc_12631_hw.m2))
row.names(acc_12631_hw) <- c("Training M1", "Test M1", "Training M2", "Test M2")
# Creating a table for combined accuracy results
# kable(acc_12631_hw, caption = "Out of Sample Performance for Self selected Model V/s Auto selected Model (HW) Store 12631", format = "latex", booktabs = T) %>% kable_styling() %>% group_rows("Self selected", 1, 2, latex_gap_space = "2em") %>% group_rows("Automatic", 3, 4, latex_gap_space = "2em")
```
```{r 12631 Holt Winter Training on all data and predicting, include=FALSE}
# Training on all data
ts_12631.ets <- ets(ts_12631, model = "AAA")
# Forecast for next two weeks
ts_12631.ets.f <- forecast(ts_12631.ets, h = 14)
# Predictions output formatting
ts_12631.ets.f.df <- as.data.frame(ts_12631.ets.f)
final_forecast_12631_hw <- cbind(days, ts_12631.ets.f.df$`Point Forecast`)
# Final Forecast
# kable(final_forecast_12631_hw, caption = "2 Week Prediction for Store 12631 using HW model", col.names = c("Day", "Forecast for Store 12631"), align = 'c')
```
```{r 12631 2 week forecast, include=FALSE}
# par(mfrow=c(2,1))
# plot(ts_12631.arima.f)
# plot(ts_12631.ets.f)
```
```{r 12631 2 week forecast comparison, include=FALSE}
forecast_12631 <- cbind(days, ts_12631.arima.f.df$`Point Forecast`, ts_12631.ets.f.df$`Point Forecast`)
colnames(forecast_12631) <- c("Day", "Arima Model Forecast", "HW Model Forecast")
forecast_12631['diff'] <- forecast_12631$`Arima Model Forecast` - forecast_12631$`HW Model Forecast`
# kable(forecast_12631, caption = "Prediction Comparison", align = 'c')
```
```{r 12631 Out of sample performance comparison, include=FALSE}
acc_12631_final <- as.data.frame(rbind(acc_12631_arima.m1, acc_12631_hw.m2))
row.names(acc_12631_final) <- c("Training M1", "Test M1", "Training M2", "Test M2")
# kable(acc_12631_final, caption = "Out of Sample Performance Comparison ARIMA v/s Holt-Winters for Store 12631", format = "latex", booktabs = T) %>% kable_styling() %>% group_rows("ARIMA", 1, 2, latex_gap_space = "2em") %>% group_rows("Holt-Winters", 3, 4, latex_gap_space = "2em")
```
```{r Store 4904, include=FALSE}
# Creating time series object
ts_4904 <- ts((data_4904$total_lettuce), frequency = 7, start = c(11, 2))
# Splitting into train and test
ts_4904_train <- subset(ts_4904, end = 85)
ts_4904_test <- subset(ts_4904, start = 86)
```
```{r 4904 Plotting, include=FALSE}
# Plot
# autoplot(ts_4904_train)
# Seasonality plot
# ggseasonplot(ts_4904_train)
# Sub series plot
# ggsubseriesplot(ts_4904_train)
```
```{r 4904 ARIMA Stationary Test, include=FALSE}
# stationary test
# Dickey Fuller Test
adf_4904 <- adf.test(log(ts_4904_train))
# Rejected the null that the sereis contain a unit root in favour of stationarity
# Phillips - Perron Test
pp_4904 <- pp.test(log(ts_4904_train))
# Rejected the null that the sereis contain a unit root in favour of stationarity
# KPSS test
kpss_4904 <- kpss.test(log(ts_4904_train))
# Fail to reject the null that the proces is STATIONARY
# Confirming using Ndiffs
ndiffs(log(ts_4904_train))
# 0 difference in needed
# Plotting the time series to see if stationary
# autoplot(log(ts_4904_train))
# Does not look stationary
# Checking how many seasonal differences are needed
nsdiffs(log(ts_4904_train))
# 1 seasonal difference in needed
# Creating a table to summarize first test results
test_results_4904 <- data.frame(Test=c("Dickey Fuller", "Phillips - Perron", "KPSS"),
"P Value" = c(adf_4904$p.value,pp_4904$p.value,kpss_4904$p.value),
"Null Hypothesis" = c("Unit Root", "Unit Root", "Stationary"),
Conclusion = c("Stationary", "Stationary", "Stationary"))
# Final test results
final_test_results_4904 <- test_results_4904
row.names(final_test_results_4904) <- c(1, 2, 3)
# Creating a table to present these results after taking 1st difference
# kable(final_test_results_4904, caption = "Stationary Test Results for Store 4904", format = "latex", booktabs = T, align = 'l')
# Taking seasonal difference
ts_4904_train.sdiff1 <- diff(log(ts_4904_train), differences = 1, lag = 7)
# Taking first difference
ts_4904_train.sdiff1.diff1 <- diff(ts_4904_train.sdiff1, differences = 1)
# Plot after seasonal diff
# autoplot(ts_4904_train.sdiff1)
# Plot after seasonal and first diff
# autoplot(ts_4904_train.sdiff1.diff1)
```
```{r 4904 ARIMA Plotting time series data, include=FALSE}
# Model 1 with only seasonal difference
# PACF AND ACF
# ggtsdisplay(ts_4904_train.sdiff1, lag.max = 49)
# Looking at the ACF and PACF model
# p = 0, q = 0 (Since it seems like white noise)
# Seasonal - ACF tails off and PACF tails off Q = 1, P = 1
# Model 2 with both seasonal and first difference
# PACF AND ACF
# ggtsdisplay(ts_4904_train.sdiff1.diff1, lag.max = 49)
# Looking at the ACF and PACF model
# p = 1, q = 1 (ACF tails off PACF also tails off)
# Seasonal - ACF cuts off at 1 and PACF tails off
# P = 0, Q = 1
```
```{r 4904 ARIMA Chosing best model, results='asis', include=FALSE}
auto.arima(ts_4904_train, trace = TRUE, ic = 'bic', lambda = 0)
# Best model is ARIMA(0,0,0)(2,1,0) (BIC - -27.06)
# Using PACF AND ACF Model we get ARIMA(0,0,0)(1,1,1)
# Using this we get a BIC -31.8 which is clearly lower than -27.06 and therefore this is
# our best model. Let's use both models on the test set to check performance
# Model 1 with only seasonal differ
ts_4904_train.arima.m1 <- Arima(ts_4904_train, order = c(0, 0, 0),
seasonal = list(order = c(1, 1, 1), period = 7), lambda = 0)
# Model 2 with seasonal and first differnce
ts_4904_train.arima.m2 <- Arima(ts_4904_train, order = c(1, 1, 1),
seasonal = list(order = c(0, 1, 1), period = 7), lambda = 0)
# Auto Arima model
ts_4904_train.arima.m3 <- Arima(ts_4904_train, order = c(0, 0, 0),
seasonal = list(order = c(2, 1, 0), period = 7), lambda = 0)
# texreg(list(ts_4904_train.arima.m1, ts_4904_train.arima.m2, ts_4904_train.arima.m3), caption = "ARIMA Two Self selected V/s automatically selected model for Store 4904",custom.model.names = c("Self (only seasonal)", "Self (both 1st and seasonal diff)", "Auto Arima"))
# Order of performance 2, 1, 3
```
```{r 4904 ARIMA Checking residual, include=FALSE}
# Model 1
# checkresiduals(ts_4904_train.arima.m1, test=FALSE)
# Model 2
# checkresiduals(ts_4904_train.arima.m2, test=FALSE)
# Model 3
# checkresiduals(ts_4904_train.arima.m3, test=FALSE)
# For both model the residual looks fine and the p-value for Ljung-Box test is way higher
# than standard 0.05 (m1 - .2764, m2 - .4851)
# Ljung-Box test
ljung_box_test_4904 <- data.frame(Model = c("ARIMA (0,0,0)(1,1,1)", "ARIMA(1,1,1)(0,1,1)", "ARIMA (0,0,0)(2,1,0)"),
"p-value" = c(0.005198, 0.1241, 0.009388))
# Print output
# kable(ljung_box_test_4904, caption = "Ljung-Box test result for the three candidate models (Store 4904)", format = "latex", booktabs = T, align = 'l')
```
```{r 4904 ARIMA accuracy test, include=FALSE}
# Accuracy test for model 1
acc_4904_arima.m1 <- accuracy(forecast(ts_4904_train.arima.m1, h = 10), ts_4904_test)
# Accuracy test for model 2
acc_4904_arima.m2 <- accuracy(forecast(ts_4904_train.arima.m2, h = 10), ts_4904_test)
# Accuracy test for best aacording to auto arima test
acc_4904_arima.m3 <- accuracy(forecast(ts_4904_train.arima.m3, h = 10), ts_4904_test)
# Creating a dataframe and renaming the rows to avoid duplicates
acc_4904_arima <- as.data.frame(rbind(acc_4904_arima.m1, acc_4904_arima.m2, acc_4904_arima.m3))
row.names(acc_4904_arima) <- c("Training M1", "Test M1", "Training M2", "Test M2", "Training M3", "Test M3")
# Creating a table for combined accuracy results
# kable(acc_4904_arima, caption = "Out of Sample Performance for The three models", format = "latex", booktabs = T) %>% kable_styling() %>% group_rows("Self selected M1", 1, 2, latex_gap_space = "2em") %>% group_rows("Self selected M2", 3, 4, latex_gap_space = "2em")%>% group_rows("Auto Arima", 5, 6, latex_gap_space = "2em")
# Model 1 outperforms all other models here However since the residual analysis
# was in favour of Model 2 We use that.
```
```{r 4904 ARIMA Training on all data and predicting, include=FALSE}
# Training on all data
ts_4904.arima.m <- Arima(ts_4904_train, order = c(1, 1, 1),
seasonal = list(order = c(0, 1, 1), period = 7), lambda = 0)
# Forecast for next two weeks
ts_4904.arima.f <- forecast(ts_4904.arima.m, h = 14)
# Final trained model output
texreg(ts_4904.arima.m,
caption = "Final Model for Store 4904 (Trained on all data)")
# Predictions output formatting
ts_4904.arima.f.df <- as.data.frame(ts_4904.arima.f)
final_forecast_4904 <- cbind(days, ts_4904.arima.f.df$`Point Forecast`)
# Final Forecast
# kable(final_forecast_4904, caption = "2 Week Prediction for Store 4904", col.names = c("Day", "Forecast for Store 4904"), align = 'c')
```
```{r 4904 Holt Winters Model, results='asis', include=FALSE}
# Looking at the range bar plots we can see additive seasonality and linear(additive)
# Range bar plots
# ts_4904_train %>% stl(s.window = "period") %>% autoplot
# Seasonality and trend has a magnitude which is less than the stochastic component.
# Seasonality ranges from +50 to -65. Indicating some days have higher sales than others
# Trend shows a continuos rise till week 15. Followed by a drastic dip in sales in week
# 20 to 21.
# trend
ts_4904.ets1 <- ets(ts_4904_train, model = "AAA")
# Automatic approach
ts_4904.ets2 <- ets(ts_4904_train, model = "ZZZ")
# According to BIC model 2 outperforms model 1
# Output texreg
# Output
# texreg(list(ts_4904.ets1, ts_4904.ets2), caption = "Holt Winters Self selected(A,A,A) V/s automatically selected model(A,N,A) for Store 4904",model.names = c("Self", "Automatic"), single.row = T)
```
```{r 4904 Holt Winters Model Out of sample accuracy, include=FALSE}
# Forecast
# Model 1
ts_4904.ets1.f <- forecast(ts_4904.ets1, h = 10)
# Model 2
ts_4904.ets2.f <- forecast(ts_4904.ets2, h = 10)
# Out of sample accuracy
# Model 1
acc_4904_hw.m1 <- accuracy(ts_4904.ets1.f, ts_4904_test)
# Model 2
acc_4904_hw.m2 <- accuracy(ts_4904.ets2.f, ts_4904_test)
# Model 2 performs better according to out of sample accuracy with most of the accuracy
# measures favouring it. Thus we go ahead with Model 2.
# Creating a dataframe and renaming the rows to avoid duplicates
acc_4904_hw <- as.data.frame(rbind(acc_4904_hw.m1, acc_4904_hw.m2))
row.names(acc_4904_hw) <- c("Training M1", "Test M1", "Training M2", "Test M2")
# Creating a table for combined accuracy results
# kable(acc_4904_hw, caption = "Out of Sample Performance for Self selected Model V/s Auto selected Model (HW) Store 4904", format = "latex", booktabs = T) %>% kable_styling() %>% group_rows("Self selected", 1, 2, latex_gap_space = "2em") %>% group_rows("Automatic", 3, 4, latex_gap_space = "2em")
```
```{r 4904 Holt Winter Training on all data and predicting, include=FALSE}
# Training on all data
ts_4904.ets <- ets(ts_4904, model = "ZZZ")
# Forecast for next two weeks
ts_4904.ets.f <- forecast(ts_4904.ets, h = 14)
# Predictions output formatting
ts_4904.ets.f.df <- as.data.frame(ts_4904.ets.f)
final_forecast_4904_hw <- cbind(days, ts_4904.ets.f.df$`Point Forecast`)
# Final Forecast
# kable(final_forecast_4904_hw, caption = "2 Week Prediction for Store 4904 using HW model", col.names = c("Day", "Forecast for Store 4904"), align = 'c')
```
```{r 4904 2 week forecast, fig.align = "center", include=FALSE}
# par(mfrow=c(2,1))
# plot(ts_4904.arima.f)
# plot(ts_4904.ets.f)
```
```{r 4904 2 week forecast comparison, include=FALSE}
forecast_4904 <- cbind(days, ts_4904.arima.f.df$`Point Forecast`, ts_4904.ets.f.df$`Point Forecast`)
colnames(forecast_4904) <- c("Day", "Arima Model Forecast", "HW Model Forecast")
forecast_4904['diff'] <- forecast_4904$`Arima Model Forecast` - forecast_4904$`HW Model Forecast`
# kable(forecast_4904, caption = "Prediction Comparison 4904", align = 'c')
```
```{r 4904 Out of sample performance comparison, include=FALSE}
acc_4904_final <- as.data.frame(rbind(acc_4904_arima.m2, acc_4904_hw.m2))
row.names(acc_4904_final) <- c("Training M1", "Test M1", "Training M2", "Test M2")
# kable(acc_4904_final, caption = "Out of Sample Performance Comparison ARIMA v/s Holt-Winters for Store 4904", booktabs = T) %>% kable_styling() %>% group_rows("ARIMA", 1, 2, latex_gap_space = "2em") %>% group_rows("Holt-Winters", 3, 4, latex_gap_space = "2em")
```
```{r Store 20974, include=FALSE}
# Creating time series object
ts_20974 <- ts((data_20974$total_lettuce), frequency = 7, start = c(12, 2))
# Splitting into train and test
# Dropping inconsistent data with start
ts_20974_train <- subset(ts_20974, end = 84, start = 7)
ts_20974_test <- subset(ts_20974, start = 85)
```
```{r 20974 Plotting, include=FALSE}
# Plot
# autoplot(ts_20974_train)
# Seasonality plot
# ggseasonplot(ts_20974_train)
# Sub series plot
# ggsubseriesplot(ts_20974_train)
```
```{r 20974 ARIMA Stationary Test, include=FALSE}
# stationary test
# Dickey Fuller Test
adf_20974 <- adf.test(log(ts_20974_train))
# Rejected the null that the sereis contain a unit root in favour of stationarity
# Phillips - Perron Test
pp_20974 <- pp.test(log(ts_20974_train))
# Rejected the null that the sereis contain a unit root in favour of stationarity
# KPSS test
kpss_20974 <- kpss.test(log(ts_20974_train))
# Fail to reject the null that the proces is STATIONARY at 1% cf level
# Confirming using Ndiffs
ndiffs(log(ts_20974_train))
# 0 difference in needed
# Confirming using NsDiffs
nsdiffs(log(ts_20974_train))
# 0 diff is needed
# Creating a table to summarize first test results
test_results_20974 <- data.frame(Test=c("Dickey Fuller", "Phillips - Perron", "KPSS"),
"P Value" = c(adf_20974$p.value,pp_20974$p.value,kpss_20974$p.value),
"Null Hypothesis" = c("Unit Root", "Unit Root", "Stationary"),
Conclusion = c("Stationary", "Stationary", "Stationary"))
# Final test results
final_test_results_20974 <- test_results_20974
row.names(final_test_results_20974) <- c(1, 2, 3)
# Creating a table to present these results
# kable(final_test_results_20974, caption = "Stationary Test Results for Store 20974", format = "latex", booktabs = T, align = 'l')
# Plotting the time series to see if stationary
# autoplot(log(ts_20974_train))
# looks like a trend
# Taking a difference
ts_20974_train.diff1 <- diff(log(ts_20974_train), differences = 1)
# Plotting the time series to see if stationary
# autoplot(ts_20974_train.diff1)
# Does look stationary
```
```{r 20974 ARIMA Plotting time series data, include=FALSE}
# Model 1 No difference
# PACF AND ACF
# ggtsdisplay(log(ts_20974_train), lag.max = 49)
# Looking at the ACF and PACF model
# p = 0, q = 0 (looks like white noise)
# Seasonal - ACF tails off and PACF tails off Q = 1, P = 1
# Model 2 - 1st diff
# PACF AND ACF
# ggtsdisplay(ts_20974_train.diff1, lag.max = 49)
# Looking at the ACF and PACF model
# p = 0, q = 1 (ACF cuts off at 1 and PACF tails off)
# Seasonal - tails off try 1- P & Q = 0 AND 1
```
```{r 20974 ARIMA Chosing best model, results = 'asis', include=FALSE}
auto.arima(ts_20974_train, trace = TRUE, ic = 'bic', lambda = 0)
# Best model is ARIMA(0,0,0)(1,0,0) (BIC - 76.54)
# Model 1
ts_20974_train.arima.m1 <- Arima(ts_20974_train, order = c(0, 0, 0),
seasonal = list(order = c(1, 0, 1), period = 7), lambda = 0)
# Model 2
ts_20974_train.arima.m2 <- Arima(ts_20974_train, order = c(0, 1, 1),
seasonal = list(order = c(1, 0, 1), period = 7), lambda = 0)
# Model 3 Auto arima
ts_20974_train.arima.m3 <- Arima(ts_20974_train, order = c(0, 0, 0),
seasonal = list(order = c(1, 0, 0), period = 7), lambda = 0)
# Order in terms of BIC - 3, 2, 1
# texreg(list(ts_20974_train.arima.m1, ts_20974_train.arima.m2, ts_20974_train.arima.m3), caption = "ARIMA Two Self selected V/s automatically selected model for Store 20974",custom.model.names = c("Self (No diff)", "Self (1st diff)", "Auto Arima"))
```
```{r 20974 ARIMA Checking residual, include=FALSE}
# Model 1
# checkresiduals(ts_20974_train.arima.m1, test=FALSE)
# Model 2
# checkresiduals(ts_20974_train.arima.m2, test=FALSE)
# Model 3
# checkresiduals(ts_20974_train.arima.m3, test=FALSE)
# p values
# model 1 - 0.7363
# model 2 - 0.8604
# model 3 - 0.7892
# Ljung-Box test
ljung_box_test_20974 <- data.frame(Model = c("ARIMA (0,0,0)(1,0,1)", "ARIMA(0,1,1)(1,0,1)", "ARIMA (0,0,0)(1,0,0)"),
"p-value" = c(0.7363, 0.8604, 0.7892))
# Print output
# kable(ljung_box_test_20974, caption = "Ljung-Box test result for the three candidate models (Store 20974)", format = "latex", booktabs = T, align = 'l')
```
```{r 20974 ARIMA accuracy test, include=FALSE}
# Accuracy for m1
acc_20974_arima.m1 <- accuracy(forecast(ts_20974_train.arima.m1, h = 10), ts_20974_test)
# Accuracy for m2
acc_20974_arima.m2 <- accuracy(forecast(ts_20974_train.arima.m2, h = 10), ts_20974_test)
# Accuracy for m3
acc_20974_arima.m3 <- accuracy(forecast(ts_20974_train.arima.m3, h = 10), ts_20974_test)
# Creating a dataframe and renaming the rows to avoid duplicates
acc_20974_arima <- as.data.frame(rbind(acc_20974_arima.m1, acc_20974_arima.m2, acc_20974_arima.m3))
row.names(acc_20974_arima) <- c("Training M1", "Test M1", "Training M2", "Test M2", "Training M3", "Test M3")
# Creating a table for combined accuracy results
# kable(acc_20974_arima, caption = "Out of Sample Performance for The three models", booktabs = T, format = "latex") %>% kable_styling() %>% group_rows("Self selected M1", 1, 2, latex_gap_space = "2em") %>% group_rows("Self selected M2", 3, 4, latex_gap_space = "2em")%>% group_rows("Auto Arima", 5, 6, latex_gap_space = "2em")
# Model three performs the best
```
```{r 20974 ARIMA Training on all data and predicting, results='asis', include=FALSE}
# Training on all data
ts_20974.arima.m <- Arima(ts_20974_train, order = c(0, 0, 0),
seasonal = list(order = c(1, 0, 0), period = 7), lambda = 0)
# Forecast for next two weeks
ts_20974.arima.f <- forecast(ts_20974.arima.m, h = 14)
# Final trained model output
# texreg(ts_20974.arima.m, caption = "Final Model for Store 20974 (Trained on all data)")
# Predictions output formatting
ts_20974.arima.f.df <- as.data.frame(ts_20974.arima.f)
final_forecast_20974 <- cbind(days, ts_20974.arima.f.df$`Point Forecast`)
# Final Forecast
# kable(final_forecast_20974, caption = "2 Week Prediction for Store 20974", col.names = c("Day", "Forecast for Store 20974"), align = 'c')
```
```{r 20974 Holt Winters Model, results='asis', include=FALSE}
# Looking at the range bar plots we can see additive seasonality and linear(additive)
# Range bar plots
# ts_20974_train %>% stl(s.window = "period") %>% autoplot
# Seasonality and trend has a magnitude which is less than the stochastic component.
# Seasonality ranges from +25 to -65. Indicating some days have higher sales than others
# Trend shows a continuos rise till week 15. Followed by a drastic dip in sales in week 16
# Further a rise till 18 and decrease after that
# trend
ts_20974.ets1 <- ets(ts_20974_train, model = "AAA")
# Automatic approach
ts_20974.ets2 <- ets(ts_20974_train, model = "ZZZ")
# Model 2 performs better in BIC terms
# Output
# texreg(list(ts_20974.ets1, ts_20974.ets2), caption = "Holt Winters Self selected(A,A,A) V/s automatically selected model(A,N,A) for Store 20974",model.names = c("Self", "Automatic"), single.row = T)
```
```{r 20974 Holt Winters Model Out of sample accuracy, include=FALSE}
# Forecast
# Model 1
ts_20974.ets1.f <- forecast(ts_20974.ets1, h = 10)
# Model 2
ts_20974.ets2.f <- forecast(ts_20974.ets2, h = 10)
# Out of sample accuracy
# Model 1
acc_20974_hw.m1 <- accuracy(ts_20974.ets1.f, ts_20974_test)
# Model 2
acc_20974_hw.m2 <- accuracy(ts_20974.ets2.f, ts_20974_test)
# Model 2 performs marginally better
# Creating a dataframe and renaming the rows to avoid duplicates
acc_20974_hw <- as.data.frame(rbind(acc_20974_hw.m1, acc_20974_hw.m2))
row.names(acc_20974_hw) <- c("Training M1", "Test M1", "Training M2", "Test M2")
# Creating a table for combined accuracy results
# kable(acc_20974_hw, caption = "Out of Sample Performance for The three models", format = "html") %>% group_rows("Self selected", 1, 2) %>% group_rows("Automatic", 3, 4) %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), position = "center", full_width = F)
```
```{r 20974 Holt Winter Training on all data and predicting, include=FALSE}
# Training on all data
ts_20974.ets <- ets(ts_20974, model = "ANA")
# Forecast for next two weeks
ts_20974.ets.f <- forecast(ts_20974.ets, h = 14)
# Predictions output formatting
ts_20974.ets.f.df <- as.data.frame(ts_20974.ets.f)
final_forecast_20974_hw <- cbind(days, ts_20974.ets.f.df$`Point Forecast`)
# Final Forecast
# kable(final_forecast_20974_hw, caption = "2 Week Prediction for Store 20974 using HW model", col.names = c("Day", "Forecast for Store 20974"), align = 'c', booktabs = T, format = "latex")
```
```{r 20974 2 week forecast, fig.align = "center", include=FALSE}
# par(mfrow=c(2,1))
# plot(ts_20974.arima.f)
# plot(ts_20974.ets.f)
```
```{r 20974 2 week forecast comparison, include=FALSE}
forecast_20974 <- cbind(days, ts_20974.arima.f.df$`Point Forecast`, ts_20974.ets.f.df$`Point Forecast`)
colnames(forecast_20974) <- c("Day", "Arima Model Forecast", "HW Model Forecast")
forecast_20974['diff'] <- forecast_20974$`Arima Model Forecast` - forecast_20974$`HW Model Forecast`
# kable(forecast_20974, caption = "Prediction Comparison 20974", align = 'c')
```
```{r 20974 Out of sample performance comparison, include=FALSE}
acc_20974_final <- as.data.frame(rbind(acc_20974_arima.m2, acc_20974_hw.m2))
row.names(acc_20974_final) <- c("Training M1", "Test M1", "Training M2", "Test M2")
# kable(acc_20974_final, caption = "Out of Sample Performance Comparison ARIMA v/s Holt-Winters for Store 20974", booktabs = T) %>% kable_styling() %>% group_rows("ARIMA", 1, 2, latex_gap_space = "2em") %>% group_rows("Holt-Winters", 3, 4, latex_gap_space = "2em")
```
```{r Store 46673, include=FALSE}
# Creating time series object
ts_46673 <- ts((data_46673$total_lettuce), frequency = 7, start = c(10, 1))
# Splitting into train and test
ts_46673_train <- subset(ts_46673, end = 93)
ts_46673_test <- subset(ts_46673, start = 94)
```
```{r 46673 Plotting, include=FALSE}
# Plot
# autoplot(ts_46673_train)
# Seasonality plot
# ggseasonplot(ts_46673_train)
# Sub series plot
# ggsubseriesplot(ts_46673_train)
```
```{r 46673 ARIMA Stationary Test, include=FALSE}
# stationary test
# Dickey Fuller Test
adf_46673 <- adf.test(log(ts_46673_train))
# Rejected the null that the sereis contain a unit root in favour of stationarity
# Phillips - Perron Test
pp_46673 <- pp.test(log(ts_46673_train))
# Rejected the null that the sereis contain a unit root in favour of stationarity
# KPSS test
kpss_46673 <- kpss.test(log(ts_46673_train))
# Fail to reject the null that the proces is STATIONARY at 1% cf level
# Confirming using Ndiffs
ndiffs(log(ts_46673_train))
# Confirming using Nsdiffs
nsdiffs(log(ts_46673_train))
# Creating a table to summarize first test results
test_results_46673<- data.frame(Test=c("Dickey Fuller", "Phillips - Perron", "KPSS"),
"P Value" = c(adf_46673$p.value,pp_46673$p.value,kpss_46673$p.value),
"Null Hypothesis" = c("Unit Root", "Unit Root", "Stationary"),
Conclusion = c("Stationary", "Stationary", "Stationary"))
# Final test results
final_test_results_46673 <- test_results_46673
row.names(final_test_results_46673) <- c(1, 2, 3)
# Creating a table to present these results after taking 1st difference
# kable(final_test_results_46673, caption = "Stationary Test Results for Store 46673", format = "latex", booktabs = T, align = 'l')
# Taking a seasonal difference
ts_46673_train.sdiff1 <- diff(log(ts_46673_train), differences = 1, lag = 7)
# Plotting the time series to see if stationary
# autoplot(ts_46673_train.sdiff1)
# Does look stationary
```
```{r 46673 ARIMA Plotting time series data, include=FALSE}
# PACF AND ACF
# ggtsdisplay(ts_46673_train.sdiff1, lag.max = 49)
# Looking at the ACF and PACF model
# p = 0, q = 0 (white noise)
# Seasonal - ACF tails off and PACF tails off at 1 Q = 1, P = 1
```
```{r 46673 ARIMA Chosing best model, include=FALSE}
auto.arima(ts_46673_train, trace = TRUE, ic = 'bic', lambda = 0)
# Best model is ARIMA(0,0,0)(0,1,1) (BIC - -25.97)
# Using PACF AND ACF Model we get ARIMA(0,0,0)(1,1,1)
ts_46673_train.arima.m1 <- Arima(ts_46673_train, order = c(0, 0, 0),
seasonal = list(order = c(1, 1, 1), period = 7), lambda = 0)
ts_46673_train.arima.m2 <- Arima(ts_46673_train, order = c(0, 0, 0),
seasonal = list(order = c(0, 1, 1), period = 7), lambda = 0)
# Model 2 performs better than Model 1 in terms of BIC
# texreg(list(ts_46673_train.arima.m1, ts_46673_train.arima.m2), caption = "ARIMA Self selected V/s automatically selected model for Store 46673",custom.model.names = c("Self", "Auto Arima"))
```
```{r 46673 ARIMA Checking residual, include=FALSE}
# Model 1
# checkresiduals(ts_46673_train.arima.m1, test=FALSE)
# Model 2
# checkresiduals(ts_46673_train.arima.m2, test=FALSE)
# p values
# model 1 - 0.145
# model 2 - 0.2038
# Ljung-Box test
ljung_box_test_46673 <- data.frame(Model = c("ARIMA(0,0,0)(1,1,1)", "ARIMA(0,0,0)(0,1,1)"),
"p-value" = c(0.145, 0.2038))
# Print output
# kable(ljung_box_test_46673, caption = "Ljung-Box test result for the three candidate models (Store 46673)", format = "latex", booktabs = T, align = 'l')
```
```{r 46673 ARIMA accuracy test, include=FALSE}
# Accuracy test for best model
acc_46673_arima.m1 <- accuracy(forecast(ts_46673_train.arima.m1, h = 10), ts_46673_test)
# Accuracy test for best aacording to auto arima test
acc_46673_arima.m2 <- accuracy(forecast(ts_46673_train.arima.m2, h = 10), ts_46673_test)
# Model 1 outperforms Model 2
# Creating a dataframe and renaming the rows to avoid duplicates
acc_46673_arima <- as.data.frame(rbind(acc_46673_arima.m1, acc_46673_arima.m2))
row.names(acc_46673_arima) <- c("Training M1", "Test M1", "Training M2", "Test M2")
# Creating a table for combined accuracy results
# kable(acc_46673_arima, caption = "Out of Sample Performance for The two models", booktabs = T, format = "latex") %>% kable_styling() %>% group_rows("Self selected", 1, 2, latex_gap_space = "2em") %>% group_rows("Auto Arima", 3, 4, latex_gap_space = "2em")
```
```{r 46673 ARIMA Training on all data and predicting, results='asis', include=FALSE}
# Training on all data
ts_46673.arima.m <- Arima(ts_46673_train, order = c(0, 0, 0),
seasonal = list(order = c(1, 1, 1), period = 7), lambda = 0)
# Forecast for next two weeks
ts_46673.arima.f <- forecast(ts_46673.arima.m, h = 14)
# Final trained model output
# texreg(ts_46673.arima.m, caption = "Final Model for Store 46673 (Trained on all data)")