Teoria

Regresja liniowa

Wiadomo pomijam:-)

Może być liczona na danych przekrojowych oraz szeregach czasowych. Oba te przypadki pokazano na przykładach.

Metoda naiwna

Prognoza jest równa ostatniej obserwacji

Metoda naiwna w wariancie sezonowym

Prognoza jest równa ostatniej obserwacji z odpowiedniego okresu (np. ten sam miesiąc z poprzedniego roku)

Metoda naiwna z dryfem

Prognoza jest równa ostatniej obserwacji plus/minus dryf; dryf jest przeciętną zmianą w zbiorze historycznych wartości, tj. \(y_{T + h} = y_T + h \frac{y_T - y_1}{T − 1}\)

Wygładzanie wykładnicze

Ważona średnia z obserwacji historycznych, tj. \(y_{T+1} = α y_{T} + α(1 − α)y_{T −1} + α(1 − α)^2 y_{T −2} + ⋯\). Gdzie \(0 ≤ α ≤ 1\). Jeżeli \(α = 1\), to \(y_{T+1} = y_{T}\), czyli liczy się tylko ostatnia obserwacja; im mniejsze \(α\) tym większe znaczenie mają obserwacje historyczne. Parametr ustala się arbitralnie. Metoda nie nadaje się do prognozowania szeregów wykazujących się trendem/sezonowością. Dla szeregów wykazujących trend można stosować metodę Holta, a trend i sezonowość Holta-Wintersa

ARIMA

Szeregi które wykazują się trendem i sezonowością określa się mianem niestacjonarnych. Szeregi niestacjonarne dają się często sprowadzić do stacjonarnych poprzez różne transformacje takie jak różnicowanie. Szereg różnic to \(y_t = y_t - y_{t-1}\).

Szereg stacjonarny jest szeregiem czysto losowym, nie zawiera żadnego komponentu systematycznego (takiego jak trend czy sezonowość). Z szeregiem stacjonarmy możemy utożsamiać biały szum, tj ciąg zmiennych losowych o jednakowym rozkładzie; wartości średniej równej zero i stałej zmienności (jednakowa wariancja). Jeżeli szereg czasowy ma być stacjonarny to jego wartości powinny osylować wokół zera a odchylenia od tegoż zera być względnie stałe:

?rnorm
set.seed(30)
y <- ts(rnorm(50))
autoplot(y) + ggtitle("White noise")

set.seed(30)
y <- ts(rnorm(50) * seq(1, 50, by=1) + seq(1, 250, by=5))
autoplot(y) + ggtitle("Not white noise")

Modele autoregresyjne: \(y_t = c + \phi_1 y_{t-1} + \phi_2 y_{t-2} + ... + \phi_p y_{t-p} + e_t\), gdzie \(e_t\) jest białym szumem. Jest to coś w rodzaju regresji ale zmiennymi objaśniającymi są opóźnione wartości zmiennej objaśnianej. Takie coś nazywa się modelem autoregresji rzędu \(p\).

Model średniej ruchomej (rzędu \(q\)): \(y_t = e_t + \phi_1 e_{t-1} + \phi_2 y_{e-2} + ... + \phi_q e_{t-q} + e_t + c\), gdzie \(e_t\) oznacza szum biały. Takie coś oznacza że \(y_t\) jest objaśniane poprzednimi błędami prognoz (które są nieobserwowalne BTW).

Połączenie różnicowania z modelami autoregresyjnymi i średniej ruchomej daje model ARIMA:

\(y_t = c + \phi_1 y_{t-1} + \phi_2 y_{t-2} + ... + \phi_p y_{t-p} + \phi_1 e_{t-1} + \phi_2 y_{e-2} + ... + \phi_q e_{t-q} + e_t\)

gdzie \(y_t'\) może być szeregiem zróżnicowanym. Takie coś nazywa się modelem \(ARIMA(p, d, q)\), gdzie \(d\) jest rzędem różnicowania.

Reasumując: tak jak w regresji liniowej czynimy pewne założenia co do modelu wg którego dane się kształtują, to podobnie jest w modelu ARIMA tylko model jest inny:

  • ARIMA(p,0,q) is an ARMA(p,q) process. The \(d\) parameter tells us how many times we need to difference the data to get a stationary trend. Thus if differencing a series once gives us an ARMA(p,q) process, we say the original series was ARIMA(p,1,q).

  • An ARMA(p,0) process is just an AR(p). Which tells us that current value today is just a weighted sum of past values. Sometimes prices are said to follow such a process.

  • An ARMA(0,q) process is MA(q) which is a moving average of past white noise. MA processes are often used to model shocks to a system (again say prices) which are gradually weakened over time.

Regresja dane przekrojowe

Dane Bank Światowy wybrane wskaźniki

Usuwamy kraje o liczbie ludności mniejszej niż 1mln

wb.indicators <- c(
  'NY.GDP.MKTP.KD',     ## GDP constant US$ 2015
  'SP.POP.TOTL',        ## Population total
  'EN.ATM.CO2E.PC',     ## CO2 emission in tons per/capita
  "EG.ELC.RNWX.ZS",     ## Electricity from renovable sources % of total
  'NY.GDP.PCAP.CD',     ## GDP per capita (current US$) 
  "NY.GDP.PCAP.PP.CD",  ## GDP per capita, PPP (current international $) 
  'EG.ELC.COAL.ZS',     ## Electricity produced from coal % of total
  'AG.PRD.LVSK.XD',     ## Livstock production index (2014--2016=100)
  "AG.CON.FERT.ZS",     ## Fertilizer consumption (kg/hectare of arable land)
  "EG.USE.PCAP.KG.OE"   ## Energy use (kg of oil equivalent per capita)
)

wb.gr.list <- read.csv("wb_groups_list.csv", sep = ';', dec = ".",   na.string="NA")
wb.gr <- wb.gr.list$code

## https://databank.worldbank.org/source/world-development-indicators
  
wb0 <- read.csv("World_Bank.csv", sep = ';', dec = ".",   na.string="NA") %>%
  filter (! code %in% wb.gr ) %>%
  pivot_wider( names_from = indicatorcode, values_from = value) %>%
  select (countryname, code, year, co2=EN.ATM.CO2E.PC, gdp=NY.GDP.PCAP.PP.CD, lstock=AG.PRD.LVSK.XD,
          fcons=AG.CON.FERT.ZS, euse=EG.USE.PCAP.KG.OE, 
          gdpt=NY.GDP.MKTP.KD,  pop=SP.POP.TOTL)
##levels(wb0$code)

Tylko kraje o liczbie ludności większej od 1mln:

wb1 <- wb0 %>% filter (pop > 1000000)

Dane dla roku 2015 dla tych krajów:

wb2015 <- wb0 %>% filter (year == 2015)
n1mln <- nrow(wb2015)

Takich krajów (obserwacji) jest 217. Regresja CO2 = a + b GDP:

model2015 <- lm (co2 ~ gdp, data=wb2015)
summary(model2015)
## 
## Call:
## lm(formula = co2 ~ gdp, data = wb2015)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.2141  -1.0334  -0.6445   0.2521  12.8007 
## 
## Coefficients:
##              Estimate Std. Error t value            Pr(>|t|)    
## (Intercept) 0.5078469  0.3268198   1.554               0.122    
## gdp         0.0002061  0.0000125  16.483 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.186 on 181 degrees of freedom
##   (34 observations deleted due to missingness)
## Multiple R-squared:  0.6002, Adjusted R-squared:  0.598 
## F-statistic: 271.7 on 1 and 181 DF,  p-value: < 0.00000000000000022
bptest(model2015)
## 
##  studentized Breusch-Pagan test
## 
## data:  model2015
## BP = 61.561, df = 1, p-value = 0.000000000000004292
durbinWatsonTest(model2015)
##  lag Autocorrelation D-W Statistic p-value
##    1     -0.05997567      2.119666   0.408
##  Alternative hypothesis: rho != 0

Obustronnie logarytmujemy (ln CO2 = a + b ln GDP):

model2015l <- lm(log(co2) ~ log(gdp), data=wb2015)
summary(model2015l)
## 
## Call:
## lm(formula = log(co2) ~ log(gdp), data = wb2015)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.3079 -0.3523 -0.0437  0.3254  1.5877 
## 
## Coefficients:
##              Estimate Std. Error t value            Pr(>|t|)    
## (Intercept) -10.29168    0.33349  -30.86 <0.0000000000000002 ***
## log(gdp)      1.18490    0.03581   33.09 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.557 on 181 degrees of freedom
##   (34 observations deleted due to missingness)
## Multiple R-squared:  0.8581, Adjusted R-squared:  0.8573 
## F-statistic:  1095 on 1 and 181 DF,  p-value: < 0.00000000000000022
bptest(model2015l)
## 
##  studentized Breusch-Pagan test
## 
## data:  model2015l
## BP = 0.078809, df = 1, p-value = 0.7789
qplot(log(gdp), log(co2), data=wb2015) + geom_smooth(method='lm')

Parametry takiego równania przy zmiennych niezależnych interpretuje się jako (częściowe jeżeli zmiennych jest więcej niż jedna) elastyczności: 1% wzrostu zmiennej niezależnej skutkuje x% wzrostu/spadku zmiennej objaśnianej. Innymi słowy wzorst nie jest stały ale tempo wzrostu jest stałe. W tym przypadku 1% wzrost GDP skutkuje 1.1% wzrostem emisji CO2.

W celu zmniejszenia zróżnicowania usuwamy małe kraje kierując się wielkością udziału w światowym PKB. Usuwamy z próby kraje, które miały 0,1% udziału w światowym GDP lub mniej (w roku 2015):

wb <- wb2015 %>%
  mutate( gdpts = gdpt/sum(gdpt, na.rm =T) * 100) %>%
  filter (gdpts >= 0.1 )

dk.n <- nrow(wb)

Takich krajów jest 66.

summary(wb)
##      countryname      code         year           co2               gdp       
##  Algeria   : 1   AGO    : 1   Min.   :2015   Min.   : 0.4588   Min.   : 3556  
##  Angola    : 1   ARE    : 1   1st Qu.:2015   1st Qu.: 3.0608   1st Qu.:12016  
##  Argentina : 1   ARG    : 1   Median :2015   Median : 5.5406   Median :26807  
##  Australia : 1   AUS    : 1   Mean   :2015   Mean   : 6.8238   Mean   :30705  
##  Austria   : 1   AUT    : 1   3rd Qu.:2015   3rd Qu.: 7.9745   3rd Qu.:46201  
##  Bangladesh: 1   BEL    : 1   Max.   :2015   Max.   :32.4706   Max.   :92968  
##  (Other)   :60   (Other):60                  NA's   :2         NA's   :1      
##      lstock           fcons               euse           gdpt               
##  Min.   : 91.57   Min.   :   4.295   Min.   :1537   Min.   :   80604080689  
##  1st Qu.: 99.35   1st Qu.:  89.204   1st Qu.:2510   1st Qu.:  185299586727  
##  Median :100.11   Median : 170.401   Median :3216   Median :  308693886543  
##  Mean   :100.10   Mean   : 304.811   Mean   :3692   Mean   : 1087170779380  
##  3rd Qu.:101.14   3rd Qu.: 266.805   3rd Qu.:4624   3rd Qu.:  836956913744  
##  Max.   :116.29   Max.   :2394.761   Max.   :7631   Max.   :18238300569000  
##                   NA's   :1          NA's   :36                             
##       pop                 gdpts        
##  Min.   :   2565708   Min.   : 0.1085  
##  1st Qu.:   9971790   1st Qu.: 0.2495  
##  Median :  31508288   Median : 0.4156  
##  Mean   :  92577427   Mean   : 1.4638  
##  3rd Qu.:  76047786   3rd Qu.: 1.1269  
##  Max.   :1379860000   Max.   :24.5572  
##

Naszą zmienną objaśnianą będzie emisja CO²; zwróćmy uwagę na duży zakres zmienności tej zmiennej. Dla ciekawości możemy wydrukować szczegółowo kto smrodzi a kto nie:

ggplot(wb, aes(x=reorder(countryname, co2), y= co2)) +
  geom_point() +
  geom_text(
    aes(label=sprintf("%.1f", co2), y= 0), size=2, color='red') +
  theme(axis.text = element_text(size = 7)) +
  coord_flip()

Wydrukuj macierz wykresów rozproszenia

## wydrukuj macierz wykresów rozproszenia (funkcja pairs)
wb4i <- wb %>% select (co2, gdp, fcons, euse)
pairs(wb4i)

Macierz korelacji

wb.corr <- wb %>% select (co2, gdp, fcons, euse) %>% na.omit() %>% cor()
kable(wb.corr)
co2 gdp fcons euse
co2 1.0000000 0.3386056 -0.0225460 0.7568814
gdp 0.3386056 1.0000000 0.2044521 0.5096831
fcons -0.0225460 0.2044521 1.0000000 -0.0055809
euse 0.7568814 0.5096831 -0.0055809 1.0000000

Regresja CO2 vs GDP

lm1 <- lm(co2 ~ gdp, data=wb)
summary(lm1)
## 
## Call:
## lm(formula = co2 ~ gdp, data = wb)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.6978 -1.5469 -0.8654  1.5096 13.3877 
## 
## Coefficients:
##               Estimate Std. Error t value       Pr(>|t|)    
## (Intercept) 1.01595586 0.93736096   1.084          0.283    
## gdp         0.00019433 0.00002556   7.603 0.000000000212 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.208 on 61 degrees of freedom
##   (3 observations deleted due to missingness)
## Multiple R-squared:  0.4866, Adjusted R-squared:  0.4782 
## F-statistic: 57.81 on 1 and 61 DF,  p-value: 0.0000000002121
qplot(gdp, co2, data=wb) + geom_smooth(method='lm')

bptest(lm1)
## 
##  studentized Breusch-Pagan test
## 
## data:  lm1
## BP = 27.282, df = 1, p-value = 0.0000001758
b.coeff.lm <- coefficients(lm1)["gdp"]

Interpretacja takiej postaci regresji jest następująca: wzrost GDP o jednostkę skutkuje wzrostem emisji CO² o 0.0001943 (w jednostkach zmiennej CO2; tonach na głowę; czyli o 0.1943347 kg)

Można powyższe pomnożyć przez tysiąc żeby nie było aż tak mało każdy wzrost GDP o 1000 USD skutkuje wzrostem emisji o 194.3346963 kg CO².

Regresja ln CO2 vs ln GDP

Model potęgowy dany jest równaniem: \(CO2 = a \cdot GDP^b\) lub w postaci równoważnej (po obustronnym zlogarytmowaniu): \(\ln(CO2) = \ln(a) + b \ln(GDP)\)

### Elasticities
lm1l <- lm(log(co2) ~ log(gdp), data=wb)
summary(lm1l)
## 
## Call:
## lm(formula = log(co2) ~ log(gdp), data = wb)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.01876 -0.33197 -0.06875  0.38150  0.90383 
## 
## Coefficients:
##             Estimate Std. Error t value             Pr(>|t|)    
## (Intercept) -7.36498    0.79554  -9.258    0.000000000000313 ***
## log(gdp)     0.89301    0.07901  11.303 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5053 on 61 degrees of freedom
##   (3 observations deleted due to missingness)
## Multiple R-squared:  0.6768, Adjusted R-squared:  0.6715 
## F-statistic: 127.8 on 1 and 61 DF,  p-value: < 0.00000000000000022
bptest(lm1l)
## 
##  studentized Breusch-Pagan test
## 
## data:  lm1l
## BP = 0.10098, df = 1, p-value = 0.7507
qplot(log(gdp), log(co2), data=wb) + geom_smooth(method='lm')

b.coeff.lml <- coefficients(lm1l)["log(gdp)"]

Interpretacja takiej postaci regresji jest następująca: 1% wzrost GDP skutkuje 0.8930136 procentowym wzrostem emisji CO² (Ekonometra inaczej s. 13).

Regresja ln CO2 vs GDP

Raczej mało sensowan w tym przypadku ale co szkodzi spróbowac: \(\ln(CO2) = \ln(a) + b \cdot GDP\)

## Semi-elasticities
lm1l2 <- lm(log(co2) ~ gdp, data=wb)
summary(lm1l2)
## 
## Call:
## lm(formula = log(co2) ~ gdp, data = wb)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.56284 -0.34638  0.03793  0.36638  1.08245 
## 
## Coefficients:
##                Estimate  Std. Error t value        Pr(>|t|)    
## (Intercept) 0.675200673 0.136635963   4.942 0.0000063630954 ***
## gdp         0.000030515 0.000003726   8.190 0.0000000000207 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6134 on 61 degrees of freedom
##   (3 observations deleted due to missingness)
## Multiple R-squared:  0.5237, Adjusted R-squared:  0.5159 
## F-statistic: 67.08 on 1 and 61 DF,  p-value: 0.0000000000207
bptest(lm1l2)
## 
##  studentized Breusch-Pagan test
## 
## data:  lm1l2
## BP = 0.28459, df = 1, p-value = 0.5937
qplot(gdp, log(co2), data=wb) + geom_smooth(method='lm')

b.coeff.slm <- coefficients(lm1l2)["gdp"]

Interpretacja takiej postaci regresji jest następująca: 1% wzrost GDP skutkuje 0,01 · 0.0000305 = 0.0030515 wzrostem emisji CO² (w jednostkach zmiennej CO2; tonach na głowę czyli 3.0514826 kg)

Regresja CO2 vs FONS (zużycie nawozów/hektar)

lm2 <- lm(co2 ~ fcons, data=wb)
summary(lm2)
## 
## Call:
## lm(formula = co2 ~ fcons, data = wb)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.430 -3.674 -1.720  1.406 25.800 
## 
## Coefficients:
##             Estimate Std. Error t value      Pr(>|t|)    
## (Intercept) 6.160713   0.890804   6.916 0.00000000301 ***
## fcons       0.002436   0.001916   1.271         0.208    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.777 on 62 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.02541,    Adjusted R-squared:  0.009688 
## F-statistic: 1.616 on 1 and 62 DF,  p-value: 0.2084
qplot(fcons, co2, data=wb) + geom_smooth(method='lm')

lm3 <- lm(co2 ~ euse, data=wb)
summary(lm3)
## 
## Call:
## lm(formula = co2 ~ euse, data = wb)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.8901 -0.8083 -0.1576  1.6217  5.2786 
## 
## Coefficients:
##              Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) 1.6222865  1.0577035   1.534     0.136    
## euse        0.0016203  0.0002644   6.128 0.0000013 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.231 on 28 degrees of freedom
##   (36 observations deleted due to missingness)
## Multiple R-squared:  0.5729, Adjusted R-squared:  0.5576 
## F-statistic: 37.55 on 1 and 28 DF,  p-value: 0.000001297
qplot(euse, co2, data=wb) + geom_smooth(method='lm')

albo (regresja wykorzystująca dane w postaci szeregów czasowych:

energy <- wb0 %>% filter (code == 'POL' & year >= 1990) %>% 
  select (year, co2, gdp, fcons, euse)

lm4 <- lm(co2 ~ gdp, data = energy)
summary(lm4)
## 
## Call:
## lm(formula = co2 ~ gdp, data = energy)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.86018 -0.35947  0.02721  0.37083  0.69556 
## 
## Coefficients:
##                Estimate  Std. Error t value             Pr(>|t|)    
## (Intercept)  8.80461445  0.18420767  47.797 < 0.0000000000000002 ***
## gdp         -0.00003643  0.00001049  -3.474              0.00175 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.451 on 27 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.3089, Adjusted R-squared:  0.2833 
## F-statistic: 12.07 on 1 and 27 DF,  p-value: 0.001748
qplot(gdp, co2, data=energy) + geom_smooth(method='lm')

lm5 <- lm(co2 ~ fcons, data=energy)
summary(lm5)
## 
## Call:
## lm(formula = co2 ~ fcons, data = energy)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.90797 -0.37515  0.08576  0.32614  0.75839 
## 
## Coefficients:
##              Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)  9.393255   0.352434  26.652 < 0.0000000000000002 ***
## fcons       -0.008355   0.002468  -3.386              0.00219 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4545 on 27 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.298,  Adjusted R-squared:  0.272 
## F-statistic: 11.46 on 1 and 27 DF,  p-value: 0.002188
lm6 <- lm(co2 ~ euse, data=energy)
summary(lm6)
## 
## Call:
## lm(formula = co2 ~ euse, data = energy)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.58375 -0.32117  0.06666  0.24584  0.63981 
## 
## Coefficients:
##               Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) -0.8453548  1.6375168  -0.516      0.61    
## euse         0.0036137  0.0006501   5.559 0.0000102 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3775 on 24 degrees of freedom
##   (5 observations deleted due to missingness)
## Multiple R-squared:  0.5628, Adjusted R-squared:  0.5446 
## F-statistic:  30.9 on 1 and 24 DF,  p-value: 0.00001018

Regresja szeregi czasowe

Zużycie energii/PC w PL

energy <- wb0 %>% filter (code == 'POL' & year >= 1990) %>% select (year, euse)

## zamiana na szereg czasowy; pierwsza obserwacja 1990
e_ts <- ts(energy$euse, start=1990)

Wykres

autoplot(e_ts) + ggtitle('Energia/PC/Polska') + 
  scale_x_continuous(breaks=seq(1990, 2030, by=2))

Dzielimy dane za zbiór uczący i testowy (training/test sets); zbiór uczący zawiera wszystkie informacje za wyjątkiem 3 ostatnich. Te 3 ostatnie obserwacje zawiera zbiór testowy:

ets.u <- window(e_ts, end=2012)
ets.t <- window(e_ts, start=2013)

autoplot(ets.u) + autolayer(ets.t)

Metody naiwne

Metody naiwne służą jako punkt odniesienia dla metod bardziej zaawansowanych. Jeżeli dokładność prognoz dla metod zaawansowanych nie różni się od dokładności uzyskanych metodami naiwnym, to oznacza to że te zaawansowane metody nic nie wnoszą.

Błądzenie losowe

hmx <- 3 ## horyzont

## rwf
f.naive  <- naive(ets.u, h=hmx)
f.naive.forecast <- forecast(f.naive, h=hmx)

e.residuals1 <- f.naive$residuals
e.fitted1 <- f.naive$fitted

##plot(f.naive)
##autoplot(f.naive)

e.acc.naive <- accuracy(f.naive.forecast, ets.t)
#e.acc.naive
autoplot(e_ts, series='empir') + autolayer(f.naive$mean, series = "f")

Naiwna z dryfem

## ?rwf
## naive with drift
f.rwf  <- rwf(ets.u, drift = T, h=hmx)
f.rwf.forecast <- forecast(f.rwf, h=hmx)

e.residuals1 <- f.rwf$residuals
e.fitted1 <- f.rwf$fitted

#plot(f.rwf)
##autoplot(f.rwf)

e.acc.rwf <- accuracy(f.rwf.forecast, ets.t)
#e.acc.rwf
autoplot(e_ts, series='empir') + autolayer(f.rwf$mean, series = "f")

Wygładzanie wykładnicze

f.es <- ses(ets.u, h=hmx)
##summary(f.es)
acc.ses <- accuracy(f.es, ets.t)

autoplot(e_ts, series='empir') + autolayer(f.es$mean, series = "f")

Metoda Holta

W dwóch wariantach (z ręcznie wybraną wartością alfa i automatycznie)

## Ustalone alpha
f.holt.1 <- holt(ets.u, h=hmx, alpha=.33)
f.holt.2 <- holt(ets.u, h=hmx)

##summary(f.holt.1)
##summary(f.holt.2)

acc.h <- accuracy(f.holt.2, ets.t)

autoplot(e_ts, series='empir') +
  autolayer(f.holt.1$mean, series='forecast1', size=1) +
  autolayer(f.holt.1$fitted, series="fitted1") +
  autolayer(f.holt.2$fitted, series="fitted2") +
  #autolayer(ets.t, series="t") +  
  autolayer(f.holt.2$mean, series='forecast2', size=1)

Arima

f.aa <- auto.arima(ets.u)
##summary(f.aa)
f.aa.forecast <- forecast(f.aa, h=hmx)

acc.arima <- accuracy(f.aa.forecast, ets.t)

autoplot(e_ts, series='empir') + 
  autolayer(f.aa$fitted, series="fitted") +
autolayer(f.aa.forecast$mean, series="forecast")

Regresja liniowa

f.lm <- tslm(ets.u ~ trend)

f.lm.forecast <- forecast(f.lm, h=hmx )

acc.lm <- accuracy(f.lm.forecast, ets.t)

autoplot(e_ts, series='empir') + 
  autolayer(f.lm$fitted, series="fitted") +
  autolayer(f.lm.forecast$mean, series="forecast")

Porównanie

A.table <- rbind( e.acc.naive,  e.acc.rwf, acc.ses, acc.h, acc.arima, acc.lm)

row.names(A.table) <- c('naive', 'naive/t', 'rwf', 'rwf/t', 
   'ses', 'ses/t', 'holt', 'holt/t', 'arima', 'arima/t', 'lm', 'lm/t')

A.table <- as.data.frame(A.table)
A.table <- A.table[order(A.table$RMSE),]
kable(A.table)
ME RMSE MAE MPE MAPE MASE ACF1 Theil’s U
lm/t 41.0092078 55.80201 41.00921 1.6105485 1.610548 0.5843194 -0.3076210 0.2824791
rwf/t -44.1735668 57.18343 47.78416 -1.7833455 1.924073 0.6808522 -0.3340916 1.0666691
arima/t -46.8273792 60.02499 50.15236 -1.8899551 2.019550 0.7145956 -0.3062132 1.1211800
holt/t -51.0984452 62.93805 51.09845 -2.0596258 2.059626 0.7280758 -0.3262953 1.1765885
naive/t -56.8130605 69.55344 56.81306 -2.2895745 2.289575 0.8095004 -0.2736388 1.3024405
ses/t -62.5745301 74.33390 62.57453 -2.5191951 2.519195 0.8915927 -0.2736388 1.3904452
arima -7.8224187 79.93331 66.32163 -0.4086869 2.634114 0.9449832 0.0036006 NA
holt -0.1372413 83.14483 67.72941 -0.0516314 2.684195 0.9650420 0.0039205 NA
ses -5.9967649 83.36405 67.15106 -0.2853171 2.663197 0.9568013 0.0017455 NA
rwf 0.0000000 85.15991 70.18287 -0.0440728 2.780738 1.0000000 -0.0543368 NA
naive -6.3197468 85.39408 70.18287 -0.2965490 2.782453 1.0000000 -0.0543368 NA
lm 0.0000000 117.63609 103.46398 -0.2222076 4.150719 1.4742058 0.7020577 NA

Emisja CO2 w PL (Bank Światowy)

carbon <- wb0 %>% filter (code == 'POL' & year >= 1990) %>% select (year, co2)

e_ts <- ts(carbon$co2, start=1990)

Wykres

autoplot(e_ts) + ggtitle('CO2/PC/Polska') + 
  scale_x_continuous(breaks=seq(1990, 2030, by=2))

Dzielimy na zbiór uczący i testowy (trzy ostatnie lata)

ets.u <- window(e_ts, end=2015)
ets.t <- window(e_ts, start=2016)

autoplot(ets.u) + autolayer(ets.t)

Metody naiwne

Błądzenie losowe

hmx <- 3 ## horyzont prognozy

f.naive  <- naive(ets.u, h=hmx)
f.naive.forecast <- forecast(f.naive, h=hmx)

##f.naive$fitted

##autoplot(f.naive)

e.acc.naive <- accuracy(f.naive.forecast, ets.t)
##e.acc.naive
autoplot(e_ts, series='empir') + autolayer(f.naive$mean, series = "f")

Błądzenie z dryfem

f.rwf  <- rwf(ets.u, drift = T, h=hmx)
f.rwf.forecast <- forecast(f.rwf, h=hmx)

e.residuals1 <- f.rwf$residuals
e.fitted1 <- f.rwf$fitted

plot(f.rwf)

autoplot(f.rwf)

e.acc.rwf <- accuracy(f.rwf.forecast, ets.t)
#e.acc.rwf

Wygładzanie wykładnicze

f.es <- ses(ets.u, h=hmx)
##summary(f.es)
acc.ses <- accuracy(f.es, ets.t)

autoplot(e_ts, series='empir') + autolayer(f.es$mean, series = "f")

Metoda Holta

f.holt.1 <- holt(ets.u, h=hmx, alpha=.33)
f.holt.2 <- holt(ets.u, h=hmx)

##summary(f.holt.1)
##summary(f.holt.2)

acc.h <- accuracy(f.holt.2, ets.t)

autoplot(e_ts, series='empir') +
  autolayer(f.holt.1$mean, series='forecast1', size=1) +
  autolayer(f.holt.1$fitted, series="fitted1") +
  autolayer(f.holt.2$fitted, series="fitted2") +
  #autolayer(ets.t, series="t") +  
  autolayer(f.holt.2$mean, series='forecast2', size=1)

Arima

f.aa <- auto.arima(ets.u)
##summary(f.aa)
f.aa.forecast <- forecast(f.aa, h=hmx)

acc.arima <- accuracy(f.aa.forecast, ets.t)

autoplot(e_ts, series='empir') + 
  autolayer(f.aa$fitted, series="fitted") +
  autolayer(f.aa.forecast$mean, series="forecast")

Regresja liniowa

f.lm <- tslm(ets.u ~ trend)

f.lm.forecast <- forecast(f.lm, h=hmx )

acc.lm <- accuracy(f.lm.forecast, ets.t)

autoplot(e_ts, series='empir') + 
  autolayer(f.lm$fitted, series="fitted") +
  autolayer(f.lm.forecast$mean, series="forecast")

Porównanie

A.table <- rbind( e.acc.naive,  e.acc.rwf, acc.ses, acc.h, acc.arima, acc.lm)

row.names(A.table) <- c('naive', 'naive/t', 'rwf', 'rwf/t', 
                        'ses', 'ses/t', 'holt', 'holt/t', 'arima', 'arima/t', 'lm', 'lm/t')

A.table <- as.data.frame(A.table)
A.table <- A.table[order(A.table$RMSE),]
kable(A.table)
ME RMSE MAE MPE MAPE MASE ACF1 Theil’s U
holt -0.0193492 0.2428023 0.1835376 -0.2620868 2.237793 0.9321091 0.0778129 NA
rwf 0.0000000 0.2467066 0.1873684 -0.0294114 2.284317 0.9515643 0.0795199 NA
arima -0.0604245 0.2497362 0.1896858 -0.7676733 2.323429 0.9633334 0.0827386 NA
ses -0.0599875 0.2497663 0.1901317 -0.7629311 2.328287 0.9655978 0.0832836 NA
naive -0.0632090 0.2546753 0.1969057 -0.8023802 2.412366 1.0000000 0.0795199 NA
lm 0.0000000 0.3538403 0.2975627 -0.1825207 3.650671 1.5111938 0.7564387 NA
naive/t 0.5130813 0.5378729 0.5130813 6.2793422 6.279342 2.6057211 -0.1697281 2.529520
arima/t 0.5130813 0.5378729 0.5130813 6.2793422 6.279342 2.6057211 -0.1697281 2.529520
ses/t 0.5130905 0.5378817 0.5130905 6.2794556 6.279456 2.6057679 -0.1697281 2.529558
holt/t 0.6036060 0.6340816 0.6036060 7.3863770 7.386377 3.0654573 -0.1172409 2.983203
rwf/t 0.6394994 0.6723561 0.6394994 7.8253186 7.825319 3.2477444 -0.1025964 3.163736
lm/t 0.6852580 0.7144489 0.6852580 8.3901891 8.390189 3.4801331 -0.1082324 3.333515

Liczba małżeństw GUS (miesięczne/BDL)

s0 <- read.csv("GUS_ludnosc.csv", sep = ';', dec = ".",   na.string="NA") 
e_ts <- ts(s0$ml, start=c(2010,1), frequency = 12)

Wykres

autoplot(e_ts)

Dzielimy na zbiór uczący i testowy (trzy ostatnie lata)

ets.u <- window(e_ts, end=c(2021,3))
ets.t <- window(e_ts, start=c(2021,4))

length(e_ts)
## [1] 143
hmx <- 8

Metody naiwne

Błądzenie losowe (wariant sezonowy)

f.naive  <- snaive(ets.u, h=hmx)
f.naive.forecast <- forecast(f.naive, h=hmx)

##f.naive$fitted

##autoplot(f.naive)

e.acc.naive <- accuracy(f.naive.forecast, ets.t)
##e.acc.naive
autoplot(e_ts, series='empir') + autolayer(f.naive$mean, series = "f")

Naiwna z dryfem

f.rwf  <- rwf(ets.u, drift = T, h=hmx)
f.rwf.forecast <- forecast(f.rwf, h=hmx)

e.residuals1 <- f.rwf$residuals
e.fitted1 <- f.rwf$fitted

plot(f.rwf)

autoplot(f.rwf)

e.acc.rwf <- accuracy(f.rwf.forecast, ets.t)
## albo
#e.acc.rwf <- accuracy(f.rwf, ets.t)
#e.acc.rwf

Wygładzanie wykładnicze

f.es <- ses(ets.u, h=hmx)
##summary(f.es)
acc.ses <- accuracy(f.es, ets.t)

autoplot(e_ts, series='empir') + autolayer(f.es$mean, series = "f") +
ggtitle('Model: SES')

Metoda Holta-Wintersa (uwzględnia sezonowość)

f.holt.1 <- hw(ets.u, h=hmx, alpha=.33)
f.holt.2 <- hw(ets.u, h=hmx)

##summary(f.holt.1)
##summary(f.holt.2)

acc.h <- accuracy(f.holt.2, ets.t)

autoplot(e_ts, series='empir') +
  autolayer(f.holt.1$mean, series='forecast1', size=1) +
  autolayer(f.holt.1$fitted, series="fitted1") +
  autolayer(f.holt.2$fitted, series="fitted2") +
  #autolayer(ets.t, series="t") +  
  autolayer(f.holt.2$mean, series='forecast2', size=1) +
  ggtitle('Model: Holt-Winters (HW)')

Arima

f.aa <- auto.arima(ets.u)
##summary(f.aa)
f.aa.forecast <- forecast(f.aa, h=hmx)

acc.arima <- accuracy(f.aa.forecast, ets.t)

autoplot(e_ts, series='empir') + 
  autolayer(f.aa$fitted, series="fitted") +
  autolayer(f.aa.forecast$mean, series="forecast") +
  ggtitle('Model: Arima (AA)')

Regresja liniowa (z sezonowością)

f.lm <- tslm(ets.u ~ trend + season)

f.lm.forecast <- forecast(f.lm, h=hmx )

acc.lm <- accuracy(f.lm.forecast, ets.t)

autoplot(e_ts, series='empir') + 
  autolayer(f.lm$fitted, series="fitted") +
  autolayer(f.lm.forecast$mean, series="forecast") +
ggtitle('Model: Linear regression (LM)')

Porównanie

A.table <- rbind( e.acc.naive,  e.acc.rwf, acc.ses, acc.h, acc.arima, acc.lm)

row.names(A.table) <- c('naive', 'naive/t', 'rwf', 'rwf/t', 
                        'ses', 'ses/t', 'holt', 'holt/t', 'arima', 'arima/t', 'lm', 'lm/t')

A.table <- as.data.frame(A.table)
A.table <- A.table[order(A.table$RMSE),]
kable(A.table)
ME RMSE MAE MPE MAPE MASE ACF1 Theil’s U
arima -0.0821922 2.045700 1.378489 -1.3937803 14.34283 0.7589711 0.0004044 NA
holt -0.0079117 2.088531 1.537541 -3.7131664 16.19963 0.8465425 0.0610292 NA
lm 0.0000000 2.138012 1.559824 -3.3762201 16.30877 0.8588108 0.3220988 NA
lm/t -0.1525899 2.239515 2.052273 -0.5388357 16.24325 1.1299443 -0.0055719 0.1931165
arima/t 0.0467939 2.360586 2.169629 -0.7214082 15.32286 1.1945583 0.0728197 0.2236221
holt/t -0.7739877 2.389823 2.099790 -5.3843061 17.17859 1.1561064 -0.0289838 0.2143462
naive -0.6861789 2.793430 1.816260 -9.0095446 18.50722 1.0000000 0.0741879 NA
naive/t 3.0375000 5.562486 4.587500 19.9540989 28.50626 2.5257945 0.5977253 0.6983175
ses -0.0253360 8.036372 6.124363 -23.0133669 55.25006 3.3719635 -0.0034137 NA
rwf 0.0000000 8.067261 6.199577 -22.6439934 55.71219 3.4133748 -0.0196408 NA
ses/t 13.6575959 16.617759 13.657596 62.7354199 62.73542 7.5196253 0.3508459 1.4135733
rwf/t 13.7641791 16.706859 13.764179 63.6005443 63.60054 7.5783081 0.3515983 1.4194806