Average cost of production
2025-11-22
library(forecast)## Warning: pakiet 'forecast' został zbudowany w wersji R 4.5.1## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(tseries)## Warning: pakiet 'tseries' został zbudowany w wersji R 4.5.1library(astsa)## Warning: pakiet 'astsa' został zbudowany w wersji R 4.5.2##
## Dołączanie pakietu: 'astsa'## Następujący obiekt został zakryty z 'package:forecast':
##
## gasMonth_Value=read.csv("C:/Users/akosc/OneDrive/Pulpit/Ania/studia/time series/Month_Value_1.csv")AC_TS<-ts(Month_Value[4], start=c(2015,1), frequency = 12)
AC_TS<- window(AC_TS, end=c(2020,4))
plot.ts(AC_TS,main = "Average cost (Jan 2015 – Apr 2020)", ylab = "Average_cost", xlab="Year",col='magenta');abline(h=mean(AC_TS),col="blue")
par(mfrow=c(2,2))
plot.ts(AC_TS,main = "Average cost (Jan 2015 – Apr 2020)", ylab = "Average_cost", xlab="Year",col='magenta');abline(h=mean(AC_TS),col="blue")
monthplot(AC_TS,col="magenta")
LM=40 # lag.max
acf(AC_TS, lag.max = LM, main="")
pacf(AC_TS, lag.max = LM, main="")
D12_AC <- diff(AC_TS, lag = 12)
par(mfrow=c(2,2))
plot.ts(D12_AC,col="magenta"); abline(h=mean(D12_AC),col="blue")
monthplot(D12_AC,col="magenta")
LM=40 # lag.max
acf(D12_AC, lag.max = LM, main="",col="magenta")
pacf(D12_AC, lag.max = LM, main="",col="magenta")
DD12_AC <- diff(D12_AC, lag = 1)
par(mfrow=c(2,2))
plot.ts(DD12_AC,col="magenta"); abline(h=mean(DD12_AC),col="blue")
monthplot(DD12_AC,col="magenta")
LM=40 # lag.max
acf(DD12_AC, lag.max = LM, main="",col="magenta")
pacf(DD12_AC, lag.max = LM, main="",col="magenta")
auto.arima(AC_TS)## Series: AC_TS
## ARIMA(0,1,1)(1,0,0)[12]
##
## Coefficients:
## ma1 sar1
## -0.7688 0.3287
## s.e. 0.0875 0.1400
##
## sigma^2 = 71332: log likelihood = -441.51
## AIC=889.02 AICc=889.43 BIC=895.45mA<-sarima(AC_TS,0,1,1,1,0,0,12,no.constant = TRUE,col = 'magenta')## initial value 5.840333
## iter 2 value 5.637467
## iter 3 value 5.603577
## iter 4 value 5.573673
## iter 5 value 5.555560
## iter 6 value 5.554274
## iter 7 value 5.553854
## iter 8 value 5.553485
## iter 9 value 5.553041
## iter 10 value 5.553038
## iter 10 value 5.553038
## iter 10 value 5.553038
## final value 5.553038
## converged
## initial value 5.612704
## iter 2 value 5.592027
## iter 3 value 5.589483
## iter 4 value 5.589186
## iter 5 value 5.589184
## iter 6 value 5.589183
## iter 6 value 5.589183
## final value 5.589183
## converged
## <><><><><><><><><><><><><><>
##
## Coefficients:
## Estimate SE t.value p.value
## ma1 -0.7688 0.0875 -8.7902 0.0000
## sar1 0.3287 0.1400 2.3472 0.0222
##
## sigma^2 estimated as 69067.68 on 61 degrees of freedom
##
## AIC = 14.11148 AICc = 14.11466 BIC = 14.21354
## 
## It is clear that residuals do satisfy assumptions.
mA$ICs## AIC AICc BIC
## 14.11148 14.11466 14.21354mA$t.table## Estimate SE t.value p.value
## ma1 -0.7688 0.0875 -8.7902 0.0000
## sar1 0.3287 0.1400 2.3472 0.0222The coefficients appear to be insignificant thus model meets our expectations.
mA_2=sarima(AC_TS,1,1,0,0,1,0,12,no.constant = TRUE,col = 'magenta')## initial value 6.046523
## iter 2 value 5.870638
## iter 3 value 5.870481
## iter 4 value 5.870444
## iter 4 value 5.870444
## final value 5.870444
## converged
## initial value 5.867894
## iter 2 value 5.867865
## iter 3 value 5.867861
## iter 3 value 5.867861
## iter 3 value 5.867861
## final value 5.867861
## converged
## <><><><><><><><><><><><><><>
##
## Coefficients:
## Estimate SE t.value p.value
## ar1 -0.5465 0.118 -4.6295 0
##
## sigma^2 estimated as 124090.7 on 50 degrees of freedom
##
## AIC = 14.65203 AICc = 14.65363 BIC = 14.72779
## 
## It is clear that residuals do not satisfy assumptions.
mA_2$ICs## AIC AICc BIC
## 14.65203 14.65363 14.72779mA_2$t.table## Estimate SE t.value p.value
## ar1 -0.5465 0.118 -4.6295 0The coefficient is significant. To summarize this model do not fit well, so for forecasting we decided to choose the best model from models compared above. :)
mA_3=sarima(AC_TS,11,1,0,0,1,0,12,col = 'magenta')## initial value 6.022301
## iter 2 value 5.916353
## iter 3 value 5.658518
## iter 4 value 5.618012
## iter 5 value 5.574839
## iter 6 value 5.546031
## iter 7 value 5.516632
## iter 8 value 5.498790
## iter 9 value 5.490830
## iter 10 value 5.487564
## iter 11 value 5.485628
## iter 12 value 5.485085
## iter 13 value 5.484731
## iter 14 value 5.484471
## iter 15 value 5.484311
## iter 16 value 5.484225
## iter 17 value 5.484191
## iter 18 value 5.484180
## iter 18 value 5.484180
## iter 18 value 5.484180
## final value 5.484180
## converged
## initial value 5.596067
## iter 2 value 5.593818
## iter 3 value 5.577463
## iter 4 value 5.571700
## iter 5 value 5.569131
## iter 6 value 5.564178
## iter 7 value 5.559193
## iter 8 value 5.555518
## iter 9 value 5.552889
## iter 10 value 5.550293
## iter 11 value 5.548181
## iter 12 value 5.545462
## iter 13 value 5.542765
## iter 14 value 5.537663
## iter 15 value 5.533433
## iter 16 value 5.529043
## iter 17 value 5.526317
## iter 18 value 5.525740
## iter 19 value 5.525009
## iter 20 value 5.524078
## iter 21 value 5.523800
## iter 22 value 5.523577
## iter 23 value 5.523494
## iter 24 value 5.523445
## iter 25 value 5.523438
## iter 26 value 5.523427
## iter 27 value 5.523426
## iter 28 value 5.523423
## iter 29 value 5.523417
## iter 30 value 5.523413
## iter 31 value 5.523411
## iter 32 value 5.523411
## iter 32 value 5.523411
## iter 32 value 5.523411
## final value 5.523411
## converged
## <><><><><><><><><><><><><><>
##
## Coefficients:
## Estimate SE t.value p.value
## ar1 -0.8857 0.1128 -7.8493 0.0000
## ar2 -0.5982 0.1517 -3.9441 0.0003
## ar3 -0.2569 0.1703 -1.5088 0.1392
## ar4 -0.0528 0.1786 -0.2956 0.7691
## ar5 -0.0888 0.1641 -0.5410 0.5915
## ar6 -0.1944 0.1582 -1.2287 0.2264
## ar7 -0.3479 0.1662 -2.0941 0.0426
## ar8 -0.0769 0.1726 -0.4455 0.6583
## ar9 0.1267 0.1760 0.7201 0.4756
## ar10 0.4830 0.1565 3.0855 0.0037
## ar11 0.5523 0.1173 4.7088 0.0000
##
## sigma^2 estimated as 55238.73 on 40 degrees of freedom
##
## AIC = 14.35529 AICc = 14.48802 BIC = 14.80983
## 
mA_3$ICs## AIC AICc BIC
## 14.35529 14.48802 14.80983mA_3$t.table## Estimate SE t.value p.value
## ar1 -0.8857 0.1128 -7.8493 0.0000
## ar2 -0.5982 0.1517 -3.9441 0.0003
## ar3 -0.2569 0.1703 -1.5088 0.1392
## ar4 -0.0528 0.1786 -0.2956 0.7691
## ar5 -0.0888 0.1641 -0.5410 0.5915
## ar6 -0.1944 0.1582 -1.2287 0.2264
## ar7 -0.3479 0.1662 -2.0941 0.0426
## ar8 -0.0769 0.1726 -0.4455 0.6583
## ar9 0.1267 0.1760 0.7201 0.4756
## ar10 0.4830 0.1565 3.0855 0.0037
## ar11 0.5523 0.1173 4.7088 0.0000This model is better than the second one but worse than the first one.
model_AC <- arima(AC_TS,order = c(0,1,1), seasonal = list(order=c(1,0,0), period=12))h <- 32 # od 2020-05 do 2022-12
future <- forecast(model_AC, h =h, level = 95)
plot(future,col="magenta")
CI_future<-data.frame(future$lower,prediction=future$mean,future$upper)
CI_future<-ts(CI_future,start=c(2020,5),end = c(2022,12),frequency = 12)
CI_future## X95. prediction X95..1
## May 2020 1394.6709 1909.764 2424.857
## Jun 2020 1379.1980 1907.874 2436.550
## Jul 2020 1427.3765 1969.295 2511.214
## Aug 2020 1565.9276 2120.773 2675.619
## Sep 2020 1398.0088 1965.487 2532.965
## Oct 2020 1245.3347 1825.170 2405.005
## Nov 2020 1231.3966 1823.331 2415.266
## Dec 2020 1183.0673 1786.859 2390.650
## Jan 2021 1344.7286 1960.149 2575.569
## Feb 2021 1215.1258 1841.959 2468.791
## Mar 2021 1146.3397 1784.381 2422.423
## Apr 2021 1276.6436 1925.700 2574.757
## May 2021 1197.0354 1907.271 2617.506
## Jun 2021 1179.0072 1906.650 2634.292
## Jul 2021 1182.1952 1926.838 2671.480
## Aug 2021 1215.3631 1976.626 2737.890
## Sep 2021 1148.0571 1925.586 2703.115
## Oct 2021 1086.0048 1879.466 2672.927
## Nov 2021 1069.7819 1878.861 2687.941
## Dec 2021 1042.4713 1866.874 2691.276
## Jan 2022 1084.3862 1923.831 2763.277
## Feb 2022 1030.7606 1884.984 2739.207
## Mar 2022 997.3090 1866.059 2734.809
## Apr 2022 1029.4707 1912.509 2795.547
## May 2022 997.8865 1906.451 2815.016
## Jun 2022 981.7175 1906.247 2830.777
## Jul 2022 972.6594 1912.883 2853.106
## Aug 2022 973.5881 1929.247 2884.907
## Sep 2022 941.6212 1912.471 2883.321
## Oct 2022 911.5055 1897.312 2883.119
## Nov 2022 896.5738 1897.114 2897.653
## Dec 2022 878.1143 1893.173 2908.232pred<-sarima.for(AC_TS,n.ahead=32,0,1,1,1,0,0,12,col='magenta')