00 R Markdown Basics

knitr::opts_chunk$set(comment = NA, warning = FALSE)

01 Dotazník

Dáta

library(readxl)
dotaznik <- read_excel("/Users/nellkabenuskova/Desktop/euba/Semester VI./ekonometrinka II./01 Základy a opakovanie/dotaznik.xlsx")

Popis dát z dotazníka

Skúmame výsledky prieskumu spokojnosti zákazníkov zábavného parku. K dispozícii máme prieskum na vzorke 500 dospelých respondentov a získali sme údaje (excelovský súbor dotaznik.xlsx) v štruktúre:

  • weekend – odpoveď yes / no, podľa toho, či navštívili park cez víkend
  • num.child – počet detí, ktoré som sprevádzal v zábavnom parku
  • distance – vzdialenosť môjho bydliska od zábavného parku
  • rides – spokojnosť (od 0 do 100) s jazdami na atrakciách
  • games – spokojnosť (od 0 do 100) s hrami ponúkanými v parku
  • wait – spokojnosť (od 0 do 100) s dĺžkou čakania na atrakcie
  • clean – spokojnosť (od 0 do 100) s čistotou v parku
  • overall – celková spokojnosť (od 0 do 100) so zábavným parkom

Analýza korelácie

Medzi celkovou spokojnosťou a spokojnosťou s atrakciami je pozitívna korelácia, hodnoty vzájomných korelácií sú takéto:

cor(dotaznik$overall, dotaznik$rides)
[1] 0.5859863

Väčšiu koreláciu, ako je medzi overall a rides, vidíme medzi overall a clean.

[1] 0.5859863
          num.child distance rides games  wait clean overall
num.child      1.00    -0.01 -0.04  0.00 -0.02 -0.01    0.32
distance      -0.01     1.00 -0.01 -0.01 -0.01  0.00    0.09
rides         -0.04    -0.01  1.00  0.46  0.31  0.79    0.59
games          0.00    -0.01  0.46  1.00  0.30  0.52    0.44
wait          -0.02    -0.01  0.31  0.30  1.00  0.37    0.57
clean         -0.01     0.00  0.79  0.52  0.37  1.00    0.64
overall        0.32     0.09  0.59  0.44  0.57  0.64    1.00

Zobrazenie korelácií pomocou teplotnej mapy (heatmap):

heatmap(cor(dotaznik[-1]))

Heatmap s upraveným vizuálom:

library(corrplot)
corrplot 0.95 loaded
corrplot(cor(dotaznik[-1]),method = "color",
         col = colorRampPalette(c("blue","white", "red"))(100))

Jednoduchý model

Veľa nám o závislosti napovedia aj hodnoty na rozptylovom grafe:

# Rozptylový graf s regresnou priamkou
rov1 <- lm(overall ~ rides, data = dotaznik)
plot(overall ~ rides, data = dotaznik, xlab = "rides", ylab = "overall")
abline(rov1, col = 4)

tie sú sústredené v okolí regresnej priamky, ktorú vyjadrujeme lineárnym regresným modelom v tvare:

\[overall_i=\beta_1+\beta_2rides_i+u_i\]

Odhadnutý model má tvar:

Zobraziť model 
summary(rov1 <- lm(overall ~ rides, data = dotaznik))

Call:
lm(formula = overall ~ rides, data = dotaznik)

Residuals:
    Min      1Q  Median      3Q     Max 
-33.597 -10.048   0.425   8.694  34.699 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -94.9622     9.0790  -10.46   <2e-16 ***
rides         1.7033     0.1055   16.14   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 12.88 on 498 degrees of freedom
Multiple R-squared:  0.3434,    Adjusted R-squared:  0.3421 
F-statistic: 260.4 on 1 and 498 DF,  p-value: < 2.2e-16

Tento model vyjadruje, že celková spokojnosť sa mení so zmenou spokojnosti s atrakciami (kladné znamienko \(\beta_2\) znamená, že s nárastom rides rastie aj overall). Pri jednotkovom náraste spokojnosti s atrakciami sa zmení celková spokojnosť priemerne o hodnotu \(\beta_2\) a celková spokojnosť sa pri nulovej spokojnosti s atrakciami priemerne rovná \(\beta_1\).

Pre konkrétneho (i-tého) zákazníka sa hodnoty odlišujú o jeho špecifiká a o ďalšie faktory nezahrnuté v modeli, ktoré v sebe zahŕňa náhodná zložka. Ak chceme podľa modelu zistiť, aká je priemerná celková spokojnosť zákazníka, ktorý bol spokojný s atrakciami na úrovni 95 bodov, dosadíme hodnotu 95 za \(rides_i\) do rovnice modelu a dostaneme odhad pre jeho overall, ten sa rovná 66.85.

Testovanie významnosti parametrov

Pre parameter \(\beta_2\) zapíšeme hypotézy: \(H_0: \beta_2 = 0, H_1: \beta_2 ≠ 0\).

Kritická hodnota Studentovho rozdelenia sa rovná 1.964739.

Absolútna hodnota t štatistiky parametra 16.1378265 je väčšia ako kritická hodnota 1.964739 preto zamietame nulovú hypotézu a parameter \(\beta_2\) je štatisticky významný.

Testovanie významnosti modelu ako celku

Hypotézy: H0: \(\beta_2\) = 0, H1: \(\beta_2\) ≠ 0.

Hodnota F štatistiky |hodnota| je väčšia ako kritická hodnota |hodnota|, preto zamietame nulovú hypotézu a model ako celok je štatisticky významný.

Komplexnejší model

Predpoklad, že na skúmaný jav pôsobí iba jeden významný faktor, môžeme v skutočných situáciách považovať za nereálny. Je logické očakávať, že pridaním ďalších faktorov vysvetlíme väčší podiel variability závislej premennej.

Predpokladajme, že celková spokojnosť závisí od všetkých premenných, ktoré máme k dispozícii. Odhadneme lineárny ekonometrický model v tvare:

\[overall_i=\beta_1+\beta_2weekend_i+\beta_3num_{c}hild_i+\beta_4distance_i+\beta_5rides_i+\beta_6wait_i+\beta_7games_i+\beta_8clean_i+u_i\]

Tento model vyjadruje, že napríklad:

  • pri jednotkovom náraste spokojnosti s atrakciami (rides) sa celková spokojnosť zmení priemerne o hodnotu beta5 pri ostatných premenných nezmenených

  • pri jednotkovom náraste spokojnosti s čakaním (wait) sa celková spokojnosť zmení priemerne o hodnotu beta6 pri ostatných premenných nezmenených

  • pri jednotkovom náraste spokojnosti s čistotou (clean) sa celková spokojnosť zmení priemerne o hodnotu beta8 pri ostatných premenných nezmenených, atď.

Testovanie významnosti parametrov

Absolútne hodnoty t štatistík parametrov |vypísať| sú väčšie ako kritická hodnota hodnota |hodnota| preto zamietame nulovú hypotézu a parametre sú štatisticky významné.

Nevýznamné sú parametre |vypísať|, čo znamená neuvažovať niektoré premenné alebo modifikovať model.

Testovanie významnosti modelu ako celku

Hypotéza má tvar: H_0: β_2=β_3=⋯=β_8=0. F štatistika sa rovná |hodnota |a keďže je väčšia ako kritická hodnota |hodnota|, z toho vyplýva, že zamietame nulovú hypotézu a model ako celok je štatisticky významný.

Koeficient determinácie R2 sa rovná |hodnota|. To znamená, že |hodnota| % variability premennej overall je vysvetlených našim lineárnym modelom.

Porovnanie modelov

Ak chceme porovnávať oba odhadnuté modely s rôznym počtom vysvetľujúcich premenných a rovnakou závislou premennou využijeme korigovaný koeficient determinácie. Ten sa tu rovná |hodnota| a je väčší ako pri prvom modeli |hodnota|.

Zároveň môžeme pridať grafické porovnanie odhadov oboch modelov s reálnymi údajmi.

Aký by bol ďalší postup??? - toť otázka

Predpoklady lineárneho modelu

  1. náhodný člen má nulovú strednú hodnotu

    E(\(u_t\)) = 0; pre všetky t

  2. náhodné členy nie sú navzájom korelované

  3. rozptyl náhodného člena je konštantný

  4. vysvetĺujúce premenné sú nestochastické

  5. vysvetľujúce premenné sú lineárne nezávislé

  6. náhodné členy majú normálne rozdelenie

02 Heteroskedasticita

Heteroskedasticita = podmienka konštantného rozptylu náhodného člena nie je splnená (3. podmienka)

Homoskedasticita

Heteroskedasticita

Príčiny existencie heteroskedasticity

Dôsledky heteroskedasticity

Testovanie heteroskedasticity

Heteroskedasticita náhodných zložiek u sa prejaví na reziduáloch, ktoré sa používajú na jej analýzu

Grafické zobrazenie: priebehy reziduálov, priebehy štvorcov reziduálov

Formálne štatistické testy: Whiteov test, Breuschov-Paganov-Godfreyov test, Godfeldov-Quandov test

Dáta

Súbor údajov o predaných domoch v Stockone (Kalifornia) v priebehu dvoch mesiacov:

  • sprice: predajná cena domu

  • age: vek domu v čase predaja

  • baths: počet kúpeľní

  • beds: počet spální

  • lgelot: dom je (1), nie je (0) postavený na väčšom pozemku

  • livarea: obytná plocha domu

  • pool: je (1), nie je (0) bazén

Načítanie dát z .xlsx súboru z webu (library(httr))

# install.packages("httr")
library(httr)
library(readxl)
GET("http://www.principlesofeconometrics.com/poe5/data/excel/stockton5.xlsx",
    write_disk(tf <- tempfile(fileext = ".xlsx")))
domy <- read_excel(tf)
rm(tf)
Prehľad dát 
domy
# A tibble: 2,610 × 7
   sprice livarea  beds baths lgelot   age  pool
    <dbl>   <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>
 1 138000      17     3   2        1    97     0
 2 105700      21     4   2.5      0    18     0
 3  22000       7     2   1        0    49     0
 4 255000      30     3   3        1    23     0
 5 203000      21     4   2        1    18     0
 6 129178      16     3   2        0     2     0
 7 140000      17     3   2        0     0     0
 8 118600      14     3   2        0     4     0
 9 125000      14     3   2        0     3     0
10 145000      17     3   2.5      0     1     0
# ℹ 2,600 more rows

Scatterplot - je predpoklad, že domy s väčším počtom spální majú väčší rozptyl ceny (heteroskedasticita)

boxplot(domy$sprice ~ domy$beds)

Ekonometrický model

Lineárny model: cena apartmánu závisí od jeho veku, počtu kúpeľní, spální, veľkosti pozemku, rozlohy a faktu, či je prítomný bazén, alebo nie

Zobraziť model so všetkými premennými 
# odhad
summary(rov <- lm(log(sprice) ~ age + baths + beds + lgelot + log(livarea) + pool, data = domy))

Call:
lm(formula = log(sprice) ~ age + baths + beds + lgelot + log(livarea) + 
    pool, data = domy)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.15484 -0.13838  0.00062  0.15193  1.96225 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)   8.9798145  0.0472413 190.084  < 2e-16 ***
age          -0.0028021  0.0002623 -10.683  < 2e-16 ***
baths         0.0175939  0.0128766   1.366 0.171949    
beds         -0.0657897  0.0092264  -7.131 1.29e-12 ***
lgelot        0.2697659  0.0217706  12.391  < 2e-16 ***
log(livarea)  1.0495496  0.0240771  43.591  < 2e-16 ***
pool          0.0839790  0.0218538   3.843 0.000125 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2364 on 2603 degrees of freedom
Multiple R-squared:  0.7098,    Adjusted R-squared:  0.7092 
F-statistic:  1061 on 6 and 2603 DF,  p-value: < 2.2e-16
srov <- summary(rov)

Počet kúpeľní v tomto modeli nie je štatisticky významný, no sám o sebe ten parameter byť nevýznamný nemusí, len medzi ostatnými premennými je taká čierna ovca.

Lineárny model bez počtu kúpeľní:

Zobraziť model bez počtu kúpeľní 
summary(rov2 <- lm(log(sprice) ~ age + beds + lgelot + log(livarea) + pool, data = domy))

Call:
lm(formula = log(sprice) ~ age + beds + lgelot + log(livarea) + 
    pool, data = domy)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.17694 -0.14011 -0.00078  0.15067  1.97358 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)   8.9568508  0.0441580 202.836  < 2e-16 ***
age          -0.0029476  0.0002397 -12.295  < 2e-16 ***
beds         -0.0629666  0.0089935  -7.001 3.21e-12 ***
lgelot        0.2699948  0.0217736  12.400  < 2e-16 ***
log(livarea)  1.0690147  0.0194135  55.066  < 2e-16 ***
pool          0.0841519  0.0218570   3.850 0.000121 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2364 on 2604 degrees of freedom
Multiple R-squared:  0.7096,    Adjusted R-squared:  0.7091 
F-statistic:  1273 on 5 and 2604 DF,  p-value: < 2.2e-16

Ponecháme pôvodný model so všetkými premennými vrátane počtu kúpeľní. Vyhodenie premennej nerieši problém heteroskedasticity (tá je stále prítomná aj po vyhodení počtu kúpeľní. Assumption - obidva modely majú podobné čísla - nezmení sa významnosť, koeficienty, R2, ani správanie reziduálov.

Vizualizácia problému heteroskedasticity

Možnosti grafického zobrazenia:

  • priebeh reziduálov v závislosti od každej nezávislej premennej

  • priebeh reziduálov v závislosti od odhadnutých hodnôt y

  • priebeh štvorcov reziduálov v závislosti od každej nezávislej premennej

  • priebeh štvorcov reziduálov v závislosti od odhadnutých hodnôt y

Závislosť štvorcov reziduálov voči premennej age a premennej beds 
plot(domy$age, resid(rov)^2, pch=16, col = "blue", main = "závislosť štvorcov reziduálov voči age")

plot(domy$beds, resid(rov)^2, pch=16, col = "green", main = "závislosť štvorcov reziduálov voci beds")

Rozpätie bodov sa mení pri každom x - máme heteroskedasticitu.

Testovanie prítomnosti heteroskedasticity

  • Whiteov test: všeobecný test, pomáha identifikovať prítomnosť heteroskedasticity

  • Breuschov-Paganov-Godfreyov test

  • Godfeldov-Quandov test

posledné dva predpokladajú rôzny tvar závislosti rozptylu

Whiteov test

Overuje, či štvorce reziduálu závisia od

  • vysvetľujúcich premenných

  • štvorcov vysvetľujúcich premenných

  • vzájomných súčtov vysvetľujúcich premenných (každá s každou) (pri menšom počte independent variables)

Chí-kvadrát rozdelenie, počet stupňov voľnosti = počet vysvetľujúcich premenných v regresnom modeli

Ak je hodnota štatistiky > kritická hodnota, zamietame H0 o heteroskedasticite.

Lineárny model, ručné výpočty Whitovej štatistiky

Prvý White test: reziduály rov závisia od vysvetľujúcich premenných a od ich štvorcov. I() = interpret this literally, spočíta ako premennú, inak “^” Rko berie ako formula operátor. Ručný výpočet = počet pozorovaní * R2

summary(white_eq <- lm(resid(rov)^2 ~ age + baths + beds + lgelot + log(livarea) + pool + I(age^2) + I(baths^2) + I(beds^2) + I(lgelot^2) + I(log(livarea)^2) + I(pool^2), data = domy))

Call:
lm(formula = resid(rov)^2 ~ age + baths + beds + lgelot + log(livarea) + 
    pool + I(age^2) + I(baths^2) + I(beds^2) + I(lgelot^2) + 
    I(log(livarea)^2) + I(pool^2), data = domy)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.2044 -0.0434 -0.0252  0.0110  3.8002 

Coefficients: (2 not defined because of singularities)
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)        8.160e-01  1.240e-01   6.579 5.72e-11 ***
age                8.037e-04  3.675e-04   2.187   0.0288 *  
baths              1.683e-02  1.613e-02   1.043   0.2969    
beds              -3.790e-02  1.933e-02  -1.961   0.0500 *  
lgelot             8.729e-03  1.088e-02   0.802   0.4224    
log(livarea)      -5.589e-01  9.572e-02  -5.839 5.91e-09 ***
pool              -1.879e-02  1.077e-02  -1.744   0.0812 .  
I(age^2)          -4.645e-06  4.587e-06  -1.013   0.3114    
I(baths^2)        -2.562e-03  3.198e-03  -0.801   0.4231    
I(beds^2)          2.981e-03  2.784e-03   1.071   0.2844    
I(lgelot^2)               NA         NA      NA       NA    
I(log(livarea)^2)  1.078e-01  1.713e-02   6.292 3.66e-10 ***
I(pool^2)                 NA         NA      NA       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1157 on 2599 degrees of freedom
Multiple R-squared:  0.0429,    Adjusted R-squared:  0.03922 
F-statistic: 11.65 on 10 and 2599 DF,  p-value: < 2.2e-16
swhite_eq <- summary(white_eq)

# ručný výpočet
Wh_stat <- 2610*0.0429; Wh_stat # (summary z White_eq); 111.969
[1] 111.969
Wh_stat <- 2610*swhite_eq$r.squared

Druhý White test: reziduály rov závisia od interakcií premenných a od ich štvorcov. Počet pozorovaní * R2

summary(white_eq2 <- lm(resid(rov)^2 ~ (age + baths + beds + lgelot + log(livarea) + pool)^2 + I(age^2) + I(baths^2) + I(beds^2) + I(lgelot^2) + I(log(livarea)^2) + I(pool^2) , data = domy))

Call:
lm(formula = resid(rov)^2 ~ (age + baths + beds + lgelot + log(livarea) + 
    pool)^2 + I(age^2) + I(baths^2) + I(beds^2) + I(lgelot^2) + 
    I(log(livarea)^2) + I(pool^2), data = domy)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.2014 -0.0423 -0.0234  0.0132  3.8022 

Coefficients: (2 not defined because of singularities)
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)          9.163e-01  1.634e-01   5.607 2.28e-08 ***
age                  3.766e-03  1.207e-03   3.120  0.00183 ** 
baths                6.455e-02  5.553e-02   1.162  0.24515    
beds                 1.823e-02  4.584e-02   0.398  0.69083    
lgelot               1.212e-01  1.041e-01   1.164  0.24450    
log(livarea)        -7.755e-01  1.474e-01  -5.263 1.54e-07 ***
pool                 2.328e-01  1.143e-01   2.036  0.04183 *  
I(age^2)            -5.348e-06  5.047e-06  -1.060  0.28937    
I(baths^2)           1.077e-03  5.287e-03   0.204  0.83861    
I(beds^2)            1.239e-02  5.747e-03   2.157  0.03111 *  
I(lgelot^2)                 NA         NA      NA       NA    
I(log(livarea)^2)    1.925e-01  3.943e-02   4.883 1.11e-06 ***
I(pool^2)                   NA         NA      NA       NA    
age:baths            1.558e-04  2.979e-04   0.523  0.60103    
age:beds             2.347e-04  2.171e-04   1.081  0.27981    
age:lgelot           8.390e-04  6.000e-04   1.398  0.16211    
age:log(livarea)    -1.486e-03  5.540e-04  -2.681  0.00738 ** 
age:pool             1.632e-03  7.443e-04   2.192  0.02844 *  
baths:beds          -1.123e-03  9.662e-03  -0.116  0.90750    
baths:lgelot        -5.812e-02  2.462e-02  -2.361  0.01830 *  
baths:log(livarea)  -2.290e-02  2.677e-02  -0.855  0.39243    
baths:pool           4.089e-02  2.989e-02   1.368  0.17136    
beds:lgelot          2.763e-02  1.929e-02   1.432  0.15227    
beds:log(livarea)   -4.549e-02  2.288e-02  -1.988  0.04690 *  
beds:pool            2.926e-02  2.108e-02   1.388  0.16522    
lgelot:log(livarea) -3.506e-02  4.663e-02  -0.752  0.45221    
lgelot:pool          3.248e-02  3.311e-02   0.981  0.32677    
log(livarea):pool   -1.649e-01  5.264e-02  -3.133  0.00175 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1148 on 2584 degrees of freedom
Multiple R-squared:  0.06311,   Adjusted R-squared:  0.05405 
F-statistic: 6.962 on 25 and 2584 DF,  p-value: < 2.2e-16
# ručný výpočet (summary z white_eq2)
Wh_stat2 <- 2610*0.06311; Wh_stat2 # 164.7171
[1] 164.7171
swhite_eq2 <- summary(white_eq2)
Wh_stat2 <- 2610*swhite_eq2$r.squared
Výpočty pomocou library(skedastic) 
# install.packages("skedastic")
library(skedastic)

# bez interakcií, df = 12
rov.white <- white(rov, interactions = FALSE); rov.white
# A tibble: 1 × 5
  statistic  p.value parameter method       alternative
      <dbl>    <dbl>     <dbl> <chr>        <chr>      
1      112. 2.44e-18        12 White's Test greater    
rov.white$statistic 
[1] 111.9699
qchisq(0.95, 12)
[1] 21.02607
# s interakciami, df = 27 (parameter)
white(rov, interactions = TRUE) 
# A tibble: 1 × 5
  statistic  p.value parameter method       alternative
      <dbl>    <dbl>     <dbl> <chr>        <chr>      
1      165. 1.04e-21        27 White's Test greater    
krit_chi27 <- qchisq(0.95, 27); krit_chi27
[1] 40.11327

Test bez interakcií: white = 111.969 (cca podobné ako cez ručný výpočet), kritická hodnota = 21.03. White > kritická hodnota, zamietame H0, heteroskedasticita je prítomná.

Test s interakciami: white = 164.72, kritická hodnota 40.11. White > kritická hodnota, zamietame H0, heteroskedasticita je prítomná.

Breuschov-Paganov-Godfreyho test (Koenkerov test)

Predpokladá rôzny tvar závislosti rozptylu (od vysvetľúcich premenných v modeli)

Chí-kvadrát rozdelenie, počet stupňov voľnosti m (počet premenných vstupujúcich do testovacej štatistiky)

Koenkerov test predpokladá rovnaký tvar závislosti s rovnakým rozdelením, ale má jednoduchší výpočet

H0: rozptyl nezávisí od vysvetľúcej premennej (homoskedasticita)

BPG test heteroskedasticity

Testuje, či je rozptyl závisí od premennej age, či je rozptyl jej funkciou

# install.packages("lmtest")
library(lmtest)
# Koenker, studentize=TRUE default
bptest(rov, ~age, data = domy)

    studentized Breusch-Pagan test

data:  rov
BP = 21.93, df = 1, p-value = 2.828e-06
qchisq(0.95, 1)
[1] 3.841459
# Breusch-Pagan
bptest(rov, ~age, studentize = FALSE, data = domy)

    Breusch-Pagan test

data:  rov
BP = 49.16, df = 1, p-value = 2.359e-12
qchisq(0.95, 1)
[1] 3.841459

Koenkerov test: BP(Koenker) = 21.93 (df = 1), kritická hodnota = 3.84. BP > kritická hodnota, zamietame H0, rozptyl je funkciou premennej, heteroskedasticita prítomná.

Breusch-Paganov test: BP = 49.16, kritická hodnota = 3.84. BP > kritická hodnota, zamietame H0.

AI aby spravila BPG test cez regresiu, nie cez funkciu

Goldfeldov-Quandtov test

Zoradí pozorovania všetkých premenných podľa vybranej premennej X (order.by), usporiadanej od najmenšej po najväčšiu hodnotu. Vynechá c prostredných pozorovaní (fraction). Odhadne parametre modelov pre obidve skupiny (s menšími a väčšími hodnotami rozptylov), a vypočíta príslušné residual sum of squares pre obe skupiny. GQ = RSS1 / RSS2

Hypotézy:

  • H0: heteroskedasticita nie je prítomná

  • H1: rozptyl je násobkom štvorca premennej (h)

Fisherovo rozdelenie, počet stupňov voľnosti = 1/2 * (n-c) - k - 1 (v tomto prípade n = 2610 (počet pozorovaní), c = 0.5 * 2160 = 522 (20% pozorovaní), k = 6 (asi počet ohraničení), potom df = 1037)

library(lmtest)
gqtest(rov, fraction = 0.2, order.by = ~age, data=domy)

    Goldfeld-Quandt test

data:  rov
GQ = 1.6832, df1 = 1037, df2 = 1037, p-value < 2.2e-16
alternative hypothesis: variance increases from segment 1 to 2
qf(0.95, 1037, 1037)
[1] 1.107609

GQ = 1.6832, kritická hodnota = 1.1076. GQ > F, zamietame H0, heteroskedasticita je prítomná, rozptyl je funkciou premennej age2.

Riešenie problému heteroskedasticity

Transformácia modelu: všetky premenné modelu vydelíme štandardnou odchýlkou. Pozorovania s väčším rozptylom dostanú menšiu váhu a naopak.

  • vážená metóda najmenších štvorcov: do funkcie lm() pridáme argument weights (1/rozptyl)

  • modely odhadu funkcie váženou metódou najmenších štvorcov podľa vykonaného testu:

    • Breusch-Pagan: weights sú prevrátené hodnoty premennej age

    • Goldfeld-Quand: weights sú prevrátené hodnoty štvorca premennej age (1/age2) (ak akceptujeme alternatívnu hypotézu testu)

  • váhy nesmú byť 0

Riešenie, ak využijeme BPQ test. Keďže váhy nesmú byť 0, je potrebné odfiltrovať tieto hodnoty a neuvažovať ich, inak sa zobrazí chybové hlásenie. Váhy sú prevrátené hodnoty premennej age.

# summary(rov_wBPQ <- lm(log(sprice) ~ age + baths + beds + lgelot + log(livarea) + pool, weights = 1/age, data = domy))
summary(rov_wBPQ <- lm(log(sprice) ~ age + baths + beds + lgelot + log(livarea) + pool, weights = 1/age, data = domy[domy$age > 0,]))

Call:
lm(formula = log(sprice) ~ age + baths + beds + lgelot + log(livarea) + 
    pool, data = domy[domy$age > 0, ], weights = 1/age)

Weighted Residuals:
     Min       1Q   Median       3Q      Max 
-0.43699 -0.03622  0.00020  0.03259  0.42187 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)   8.8254169  0.0449551 196.316  < 2e-16 ***
age          -0.0041994  0.0003304 -12.712  < 2e-16 ***
baths         0.0909612  0.0131886   6.897 6.68e-12 ***
beds         -0.1171450  0.0091612 -12.787  < 2e-16 ***
lgelot        0.2111857  0.0210699  10.023  < 2e-16 ***
log(livarea)  1.1284603  0.0237918  47.431  < 2e-16 ***
pool          0.0270659  0.0215706   1.255     0.21    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.05953 on 2536 degrees of freedom
Multiple R-squared:  0.7773,    Adjusted R-squared:  0.7768 
F-statistic:  1475 on 6 and 2536 DF,  p-value: < 2.2e-16

Riešenie, ak využijeme GQ test - tiež filtrujeme nulu nadobúdajúce hodnoty premennej age. Váhy sú prevrátené hodnoty štvorca premennej.

# summary(rov_wGQ <- lm(log(sprice) ~ age + baths + beds + lgelot + log(livarea) + pool, weights = 1/age^2, data = domy))
summary(rov_wGQ <- lm(log(sprice) ~ age + baths + beds + lgelot + log(livarea) + pool, weights = 1/age^2, data = domy[domy$age > 0,]))

Call:
lm(formula = log(sprice) ~ age + baths + beds + lgelot + log(livarea) + 
    pool, data = domy[domy$age > 0, ], weights = 1/age^2)

Weighted Residuals:
     Min       1Q   Median       3Q      Max 
-0.38238 -0.00704  0.00244  0.00877  0.30968 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)   8.7693445  0.0359147 244.171  < 2e-16 ***
age          -0.0079088  0.0005949 -13.293  < 2e-16 ***
baths         0.2123953  0.0110608  19.203  < 2e-16 ***
beds         -0.1603081  0.0082846 -19.350  < 2e-16 ***
lgelot        0.1098063  0.0174078   6.308 3.33e-10 ***
log(livarea)  1.1225561  0.0199958  56.140  < 2e-16 ***
pool          0.0268165  0.0263771   1.017    0.309    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.02399 on 2536 degrees of freedom
Multiple R-squared:  0.8636,    Adjusted R-squared:  0.8633 
F-statistic:  2677 on 6 and 2536 DF,  p-value: < 2.2e-16

Baths je už štatisticky významný parameter, pool je problém (nie je štatisticky významný), parameter počtu spální beds je záporný, riešenie - kvadratická funkcia.

V korelačnej matici vidíme silné korelácie medzi log(livarea) a baths / beds. Podľa džepeta ide o multikolinearitu - keď je väčší dom, logicky v ňom bude aj viac spální a kúpeľní (teda neviem odhadnúť ich samostatný vplyv).

cor(model.matrix(rov)) 
             (Intercept)          age      baths       beds      lgelot
(Intercept)            1           NA         NA         NA          NA
age                   NA  1.000000000 -0.4981936 -0.3589394 0.003568086
baths                 NA -0.498193585  1.0000000  0.6472284 0.146827032
beds                  NA -0.358939430  0.6472284  1.0000000 0.071604302
lgelot                NA  0.003568086  0.1468270  0.0716043 1.000000000
log(livarea)          NA -0.292157132  0.7656310  0.6559962 0.228554334
pool                  NA -0.054155615  0.1651449  0.1415390 0.130282564
             log(livarea)        pool
(Intercept)            NA          NA
age            -0.2921571 -0.05415561
baths           0.7656310  0.16514486
beds            0.6559962  0.14153898
lgelot          0.2285543  0.13028256
log(livarea)    1.0000000  0.21363241
pool            0.2136324  1.00000000

Lineárny model nemusí byť vhodný. Finálny model by mohol využívať kvadratickú beds2. Váhy sú podľa Koenkerovho testu weights = 1/age. Úpravy výsledného modelu:

  • beds + I(beds^2) : nelineárny efekt počtu spální, efekt spální sa mení podľa toho, koľko ich už je (kvadratická funkcia)

  • interakcia lgelot:pool: efekt bazéna závisí od toho, či má dom veľký pozemok

summary(rov_wBPQ_final <- lm(log(sprice) ~ age + baths + beds + I(beds^2) + lgelot + log(livarea) + lgelot:pool, weights = 1/age, data=domy[domy$age>0,]))

Call:
lm(formula = log(sprice) ~ age + baths + beds + I(beds^2) + lgelot + 
    log(livarea) + lgelot:pool, data = domy[domy$age > 0, ], 
    weights = 1/age)

Weighted Residuals:
     Min       1Q   Median       3Q      Max 
-0.43178 -0.03576  0.00010  0.03258  0.42364 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)   8.9831974  0.0924136  97.206  < 2e-16 ***
age          -0.0043888  0.0003422 -12.825  < 2e-16 ***
baths         0.0855283  0.0133267   6.418 1.64e-10 ***
beds         -0.2060610  0.0461445  -4.466 8.34e-06 ***
I(beds^2)     0.0128691  0.0065600   1.962   0.0499 *  
lgelot        0.1943508  0.0226115   8.595  < 2e-16 ***
log(livarea)  1.1304709  0.0237270  47.645  < 2e-16 ***
lgelot:pool   0.1174009  0.0516886   2.271   0.0232 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.05945 on 2535 degrees of freedom
Multiple R-squared:  0.778, Adjusted R-squared:  0.7773 
F-statistic:  1269 on 7 and 2535 DF,  p-value: < 2.2e-16

Autoregresne podmienená heteroskedasticita

03 Multikolinearita

Ak je porušená podmienka lineárnej nezávislosti vysvetľujúcich premenných, gratulujem, máme prítomnú multikolinearitu.

Dôsledky prítomnosti multikolinearity

Dôsledk prítomnosti multikolinearity vysokého stupňa:

  • odhadnuté rozptyly parametrov (aj štandardné odchýlky) sú veľké, aj napriek neskreslenosti a efektívnosti odhadu parametrov

  • nie je možné úplne oddeliť vplyvy jednotlivých vysvetľujúcich premenných na závislú premennú (tým pádom ani ich správne interpretovať)

  • vysoké hodnoty rozptylov parametrov => nízke hodnoty t-štatistík, teda prijatie hypotézy o nevýznamnosti parametra

  • napriek nízkym hodnotám t-štatistík môže byť koeficient determinácie R2 veľmi vysoký

  • odhady parametrov a ich roptylov sú citlivé na malé zmeny v údajoch

Testovanie a indikátory prítomnosti multikolinearity

Multikolinearita vysokého stupňa môže byť prítomná, ak:

  • koeficient determinácie R2 je vysoký, ale niektoré t-štatistiky ukazujú nevýznamnosť parametrov

  • korelačný koeficient r > 0.9 alebo r > R2

  • pri regresii vysvetľujúcej premennej (vysvetlený ostaatnými premennými) vysoký koeficient R2

  • variance inflation factor (VIF) > 10, vypočítaný z koeficientov Rj2 (pomocné regresie)

  • Farrarov-Glauberov test

Riešenie multikolinearity

  • kombinácia prierezových údajov a časových radov

  • pri silnej závislosti dvoch vysvetľujúcich premenných jednu z nich vynecháme, a odhadujeme parametre redukovaného modelu (pozor na chyby špecifikácie)

  • transformácia premenných - diferencovanie po sebe nasledujúcich hodnôt časového radu a ich využitie v analýze (nevhodné pre prierezové údaje)

  • získanie dodatočných údajov, iný výber pozorovaní, zväčšenie rozsahu výberu (vysoký stupeň MK môže byť problémom výberu)

  • metóda skresleného odhadu parametrov: ponúka skreslený odhad, ale s menším rozptylom (napr. hrebeňová (ridge) regresia)

Príklad 1

údaje o výdavkoch na výskum a vývoj, sales a profits

Ekonometrický model

Porovnanie dvoch modelov

  • výdavky na R&D závisiace od predaja a zisku

  • výdavky na R&D závisiace od zisku

summary(rov<-lm(rd ~ sales + profits, data = ekon1))

Call:
lm(formula = rd ~ sales + profits, data = ekon1)

Residuals:
    Min      1Q  Median      3Q     Max 
-7719.4  -653.6    45.2   491.0  6590.7 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)   
(Intercept) -502.72583 1208.29714  -0.416  0.68326   
sales          0.05852    0.01491   3.925  0.00135 **
profits       -0.15779    0.17545  -0.899  0.38266   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3240 on 15 degrees of freedom
Multiple R-squared:  0.6277,    Adjusted R-squared:  0.5781 
F-statistic: 12.65 on 2 and 15 DF,  p-value: 0.0006046
summary(rov2<-lm(rd ~ profits, data = ekon1))

Call:
lm(formula = rd ~ profits, data = ekon1)

Residuals:
    Min      1Q  Median      3Q     Max 
-7282.2 -1932.9  -936.5   547.6 14170.1 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept)  796.7743  1601.9235   0.497   0.6257  
profits        0.3619     0.1587   2.281   0.0366 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4467 on 16 degrees of freedom
Multiple R-squared:  0.2454,    Adjusted R-squared:  0.1982 
F-statistic: 5.203 on 1 and 16 DF,  p-value: 0.03659

V prvom prípade rast profits zníži výdavky na R&D, v druhom prípade rast zisku výdavky na R&D naopak zvýši. Prvý regresný model rieši o koľko sa zmenia výdavky na R&D, ak zvýšim profits a sales nechám rovnaké (efekt zisku pri konštantných sales), no druhý model rieši celkový efekt zisku. Podozrenie na multikolinearitu vychádza z nemožnosti oddelenia vplyvov premenných, citlivosti na špecifikáciu a z nestabilných koeficientov.

Koeficient determinácie (1. rovnica): R2 = 0.6277

Koeficient determinácie (2. rovnica): R2 = 0.2454

Korelačná matica

Lineárna závislosť sa meria pomocou korelačných koeficientov. Absolútne hodnoty nad 0.8, prípadne nad 0.9 sú považované za vysokú lineárnu závislosť. Ak nemáme v dátach zbytočné stĺpce, použijeme cor(model.matrix(rov))

cor(model.matrix(rov)[,-1])
            sales   profits
sales   1.0000000 0.7547188
profits 0.7547188 1.0000000

Korelačný koeficient má hodnotu r = 0.7547, nie je väčší ako 0.9.

Pomocné regresie medzi vysvetľujúcimi premennými

Vysvetľujúca premenná sales vysvetlená druhou vysvetľujúcou premennou profits, koeficient determinácie porovnávame s koeficientom v prvom modeli

summary(rov.pom <- lm(sales ~ profits, data = ekon1)) 

Call:
lm(formula = sales ~ profits, data = ekon1)

Residuals:
   Min     1Q Median     3Q    Max 
-64826 -23154 -12442   6560 186721 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 22204.916  19483.927   1.140 0.271205    
profits         8.881      1.930   4.602 0.000295 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 54330 on 16 degrees of freedom
Multiple R-squared:  0.5696,    Adjusted R-squared:  0.5427 
F-statistic: 21.17 on 1 and 16 DF,  p-value: 0.0002948

R2 prvého modelu 0.6277, R2 pomocnej regresie 0.5696, nie je vyšší ako v pôvodom modeli.

Faktor inflácie rozptylu (VIF)

Vzťah na výpočet: VIF = 1 / (1 - Rj2)

Ak je VIF > 10, naznačuje to problém so silnou multikolinearitou

# install.packages("car")
library("car")
Loading required package: carData
vif(rov) 
   sales  profits 
2.323423 2.323423 
(vif.value <- 1 / (1 - summary(rov.pom)$r.squared))
[1] 2.323423

VIF = 2.3234 < 10.

Závery vyššie uvedených testovacích metód síce naznačujú, že nie je extrémna multikolinearita, ale neznamená to, že multikolinearita nie je prítomná.

Testovanie pomocou library(mctest)

Balíček na testovanie multikolinearity library(mctest) robí viac testov naraz.

omcdiag: testuje MK na úrovni celého modelu (obsahuje viac testov), 1 = áno, 0 = nie. Ak viacero testov naznačí prítomnosť MK (3 - 4), môže to byť problém.

imcdiag: individuálne testy premenných (napr. VIF)

# install.packages("mctest")
library(mctest)
omcdiag(rov)

Call:
omcdiag(mod = rov)


Overall Multicollinearity Diagnostics

                       MC Results detection
Determinant |X'X|:         0.4304         0
Farrar Chi-Square:        13.0671         1
Red Indicator:             0.7547         1
Sum of Lambda Inverse:     4.6468         0
Theil's Method:            0.5115         1
Condition Number:          4.9606         0

1 --> COLLINEARITY is detected by the test 
0 --> COLLINEARITY is not detected by the test
imcdiag(rov)

Call:
imcdiag(mod = rov)


All Individual Multicollinearity Diagnostics Result

           VIF    TOL      Wi  Fi Leamer   CVIF Klein   IND1 IND2
sales   2.3234 0.4304 21.1748 Inf  0.656 5.8855     0 0.0269    1
profits 2.3234 0.4304 21.1748 Inf  0.656 5.8855     0 0.0269    1

1 --> COLLINEARITY is detected by the test 
0 --> COLLINEARITY is not detected by the test

profits , coefficient(s) are non-significant may be due to multicollinearity

R-square of y on all x: 0.6277 

* use method argument to check which regressors may be the reason of collinearity
===================================

Pri overall testovaní 3 testy indikujú problém, 3 nie. Individuálne testovanie naznačilo problém - nie je to extrémny prípad, ale dostatočný na to, aby koeficienty boli nestabilné a interpretácia bola problém.

Príklad 2

Ekonometrický model

Summary naznačuje významnosť parametrov podľa t-štatistík, koeficient determinácie R2 = 0.951.

library(readxl)
ekon2 <- read_excel("/Users/nellkabenuskova/Desktop/euba/Semester VI./ekonometrinka II./03 Multikolinearita/ekon1.xlsx", sheet = "Kap11Tab43")
summary(rov3 <- lm(cost ~ age + miles, data = ekon2))

Call:
lm(formula = cost ~ age + miles, data = ekon2)

Residuals:
    Min      1Q  Median      3Q     Max 
-562.30 -239.24   67.94  175.13  588.03 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   26.189    114.201   0.229    0.819    
age           28.016      2.776  10.094 4.91e-14 ***
miles       -154.635     20.688  -7.475 6.99e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 259.6 on 54 degrees of freedom
Multiple R-squared:  0.951, Adjusted R-squared:  0.9492 
F-statistic: 523.8 on 2 and 54 DF,  p-value: < 2.2e-16

Korelačná matica

Korelačný koeficient r = 0.9964646 > 0.9, naznačuje vysokú MK.

cor(model.matrix(rov3)[,-1])
            age     miles
age   1.0000000 0.9964646
miles 0.9964646 1.0000000

Pomocné regresie medzi vysvetľujúcimi premennými

Koeficient determinácie R2 = 0.9929 > 0.951. Má vyššiu hodnotu ako v pôvodnom modeli - MK?

summary(rov.pom3 <- lm(miles ~ age, data = ekon2))

Call:
lm(formula = miles ~ age, data = ekon2)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.0856 -1.2822  0.3409  1.2393  2.5533 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.21719    0.48028   8.781 4.75e-12 ***
age          0.13369    0.00152  87.961  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.692 on 55 degrees of freedom
Multiple R-squared:  0.9929,    Adjusted R-squared:  0.9928 
F-statistic:  7737 on 1 and 55 DF,  p-value: < 2.2e-16

Testovanie pomocou library(mctest) a VIF

Faktor inflácie rozptylu VIF = 141.6764 > 10. Pri overall testoch má každý z nich MC detected (1), rovnako tak v individuálnych testoch. Multikolinearita by teda mala byť fest prítomná a fest problém.

# install.packages("car")
library(car)
vif(rov3)
     age    miles 
141.6764 141.6764 
library(mctest)
omcdiag(rov3) 

Call:
omcdiag(mod = rov3)


Overall Multicollinearity Diagnostics

                       MC Results detection
Determinant |X'X|:         0.0071         1
Farrar Chi-Square:       269.9682         1
Red Indicator:             0.9965         1
Sum of Lambda Inverse:   283.3528         1
Theil's Method:            1.0349         1
Condition Number:         63.7616         1

1 --> COLLINEARITY is detected by the test 
0 --> COLLINEARITY is not detected by the test
imcdiag(rov3)

Call:
imcdiag(mod = rov3)


All Individual Multicollinearity Diagnostics Result

           VIF    TOL       Wi  Fi Leamer   CVIF Klein  IND1 IND2
age   141.6764 0.0071 7737.202 Inf  0.084 -9.153     1 1e-04    1
miles 141.6764 0.0071 7737.202 Inf  0.084 -9.153     1 1e-04    1

1 --> COLLINEARITY is detected by the test 
0 --> COLLINEARITY is not detected by the test

* all coefficients have significant t-ratios

R-square of y on all x: 0.951 

* use method argument to check which regressors may be the reason of collinearity
===================================

Príklad 3

Údaje o spotrebe (GCQ) a disponibilnom príjme (GYDQ), spotrebná funkcia

library(readxl)
mydata <- read_excel("/Users/nellkabenuskova/Desktop/euba/Semester VI./ekonometrinka II./03 Multikolinearita/ekon1.xlsx", sheet = "Kap8Tab30")

Ekonometrický model

Spotrebná funkcia. Koeficient determinácie R2 = 0.9971. Prítomná je autokorelácia prvého rádu (Ljung-Box test), riešením bude oneskorenie premenných

summary(rov.spot <-lm(GCQ ~ GYDQ, data = mydata))

Call:
lm(formula = GCQ ~ GYDQ, data = mydata)

Residuals:
    Min      1Q  Median      3Q     Max 
-98.310 -26.226   3.062  31.122  78.788 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -59.7206    11.1218   -5.37 3.04e-07 ***
GYDQ          0.9340     0.0042  222.36  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 40.47 on 146 degrees of freedom
Multiple R-squared:  0.9971,    Adjusted R-squared:  0.997 
F-statistic: 4.945e+04 on 1 and 146 DF,  p-value: < 2.2e-16
# AR(1)
Box.test(resid(rov.spot), lag = 1, type = "Lj")

    Box-Ljung test

data:  resid(rov.spot)
X-squared = 106.93, df = 1, p-value < 2.2e-16

Transformácia modelu na time series. Stále prítomná autokorelácia prvého rádu. R2 = 0.9971

# install.packages("dynlm")
library(dynlm)
mydata <- ts(mydata, start = c(1958, 4), frequency = 4)
summary(rov.spot2 <- dynlm(GCQ ~ GYDQ + L(GYDQ, 1), data = mydata))

Time series regression with "ts" data:
Start = 1959(1), End = 1995(3)

Call:
dynlm(formula = GCQ ~ GYDQ + L(GYDQ, 1), data = mydata)

Residuals:
     Min       1Q   Median       3Q      Max 
-102.168  -24.010    3.737   31.787   78.558 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -57.7620    11.2126  -5.152 8.36e-07 ***
GYDQ          0.6460     0.1327   4.868 2.93e-06 ***
L(GYDQ, 1)    0.2892     0.1331   2.173   0.0314 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 40.05 on 144 degrees of freedom
Multiple R-squared:  0.9971,    Adjusted R-squared:  0.9971 
F-statistic: 2.482e+04 on 2 and 144 DF,  p-value: < 2.2e-16
# AR (1)
Box.test(resid(rov.spot2), lag = 1, type = "Lj")

    Box-Ljung test

data:  resid(rov.spot2)
X-squared = 121.79, df = 1, p-value < 2.2e-16

Korelačná matica

Z korelačnej matice vidíme vysokú koreláciu medzi premennými r = 0.9995.

cor(model.matrix(rov.spot2)[,-1])
                GYDQ L(GYDQ, 1)
GYDQ       1.0000000  0.9995005
L(GYDQ, 1) 0.9995005  1.0000000

Transformácia modelu (riešiaca multikolinearitu)

Prvá transformácia v modeli uvažuje aktuálny príjem a zmenu oproti predošlému obdobiu (d(GYDQ,1)), no evidentne nemôže ostať takto, lebo obsahuje na seba prepojené premenné (príjem a jeho zmenu), čo podľa džepeta vedie k multikolinearite. Korelácia r = 0.10617.

# prvý pokus
summary(rov.spot2t <- dynlm(GCQ ~ GYDQ + d(GYDQ, 1), data = mydata)) 

Time series regression with "ts" data:
Start = 1959(1), End = 1995(3)

Call:
dynlm(formula = GCQ ~ GYDQ + d(GYDQ, 1), data = mydata)

Residuals:
     Min       1Q   Median       3Q      Max 
-102.168  -24.010    3.737   31.787   78.558 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -57.762012  11.212567  -5.152 8.36e-07 ***
GYDQ          0.935274   0.004218 221.749  < 2e-16 ***
d(GYDQ, 1)   -0.289229   0.133088  -2.173   0.0314 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 40.05 on 144 degrees of freedom
Multiple R-squared:  0.9971,    Adjusted R-squared:  0.9971 
F-statistic: 2.482e+04 on 2 and 144 DF,  p-value: < 2.2e-16
cor(model.matrix(rov.spot2t)[,-1])
                GYDQ d(GYDQ, 1)
GYDQ       1.0000000  0.1061714
d(GYDQ, 1) 0.1061714  1.0000000

V druhej transformácii je vysvetľujúcou premennou zmena oproti minulému obdobiu a príjem z minulého obdobia (lagged one period). Korelácia r = 0.0746943

# druhý pokus
summary(rov.spot2t2 <- dynlm(GCQ ~ d(GYDQ,1) + L(GYDQ,1), data = mydata))

Time series regression with "ts" data:
Start = 1959(1), End = 1995(3)

Call:
dynlm(formula = GCQ ~ d(GYDQ, 1) + L(GYDQ, 1), data = mydata)

Residuals:
     Min       1Q   Median       3Q      Max 
-102.168  -24.010    3.737   31.787   78.558 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -57.762012  11.212567  -5.152 8.36e-07 ***
d(GYDQ, 1)    0.646044   0.132707   4.868 2.93e-06 ***
L(GYDQ, 1)    0.935274   0.004218 221.749  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 40.05 on 144 degrees of freedom
Multiple R-squared:  0.9971,    Adjusted R-squared:  0.9971 
F-statistic: 2.482e+04 on 2 and 144 DF,  p-value: < 2.2e-16
cor(model.matrix(rov.spot2t2)[,-1]) #r = 0.0746943
           d(GYDQ, 1) L(GYDQ, 1)
d(GYDQ, 1) 1.00000000 0.07469434
L(GYDQ, 1) 0.07469434 1.00000000

Tretia transformácia. Autokorelácia 1: rádu už nie je prítomná, multikolinearita stále ostáva.

summary(rov.spot3 <- dynlm(GCQ ~ GYDQ + d(GYDQ, 1) + L(GCQ, 1), data = mydata))

Time series regression with "ts" data:
Start = 1959(1), End = 1995(3)

Call:
dynlm(formula = GCQ ~ GYDQ + d(GYDQ, 1) + L(GCQ, 1), data = mydata)

Residuals:
    Min      1Q  Median      3Q     Max 
-52.111  -7.331   0.438   9.065  34.452 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.58332    4.43985   0.582   0.5616    
GYDQ         0.04709    0.02831   1.663   0.0985 .  
d(GYDQ, 1)   0.21386    0.05014   4.265 3.61e-05 ***
L(GCQ, 1)    0.95263    0.03032  31.415  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 14.3 on 143 degrees of freedom
Multiple R-squared:  0.9996,    Adjusted R-squared:  0.9996 
F-statistic: 1.301e+05 on 3 and 143 DF,  p-value: < 2.2e-16
Box.test(resid(rov.spot3), lag = 1, type = "Lj") #autokorelácia 1.radu nie je

    Box-Ljung test

data:  resid(rov.spot3)
X-squared = 0.12343, df = 1, p-value = 0.7253
cor(model.matrix(rov.spot3)[,-1]) #je multikolinearita
                GYDQ d(GYDQ, 1)  L(GCQ, 1)
GYDQ       1.0000000 0.10617137 0.99843627
d(GYDQ, 1) 0.1061714 1.00000000 0.08825099
L(GCQ, 1)  0.9984363 0.08825099 1.00000000

Štvrtá transformácia.

#transformácia na DU - poiznámky do tabletu od ADI - ja som ho mala vybitý 
summary(rov.spot3t2 <- dynlm(GCQ ~ d(GYDQ, 1) + L(GYDQ, 1) + L(GCQ, 1), data = mydata))

Time series regression with "ts" data:
Start = 1959(1), End = 1995(3)

Call:
dynlm(formula = GCQ ~ d(GYDQ, 1) + L(GYDQ, 1) + L(GCQ, 1), data = mydata)

Residuals:
    Min      1Q  Median      3Q     Max 
-52.111  -7.331   0.438   9.065  34.452 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.58332    4.43985   0.582   0.5616    
d(GYDQ, 1)   0.26095    0.04894   5.332  3.7e-07 ***
L(GYDQ, 1)   0.04709    0.02831   1.663   0.0985 .  
L(GCQ, 1)    0.95263    0.03032  31.415  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 14.3 on 143 degrees of freedom
Multiple R-squared:  0.9996,    Adjusted R-squared:  0.9996 
F-statistic: 1.301e+05 on 3 and 143 DF,  p-value: < 2.2e-16
#to moze aj pokaziť!
Box.test(resid(rov.spot3t2), lag = 1, type = "Lj")

    Box-Ljung test

data:  resid(rov.spot3t2)
X-squared = 0.12343, df = 1, p-value = 0.7253
cor(model.matrix(rov.spot3t2)[,-1])
           d(GYDQ, 1) L(GYDQ, 1)  L(GCQ, 1)
d(GYDQ, 1) 1.00000000 0.07469434 0.08825099
L(GYDQ, 1) 0.07469434 1.00000000 0.99850184
L(GCQ, 1)  0.08825099 0.99850184 1.00000000

Doplniť riešenie ridge regression a farrar glauber test

04 Časové rady, biely šum, náhodná prechádzka, stacionarita

Biely šum

Časový rad, ktorého pozorovania sú generované náhodným výberom (stochastický proces), majú pravdepodobnostné rozdelenie s priemerom E(yt) = β0 a s rozptylom var(yt) = σ2.

Stochastický proces: séria náhodných premenných usporiadaných v čase (zápis), teda sú identicky a nezávisle rozdelené (+ nejaký ďalší blábol, s ktorým sa dostal k náhodnej zložke, model a tak)

Biely šum: klasická náhodná zložka v terminológii časových radov, ak platí (model), trend je konštanta a premenná okolo nej fluktuuje (biely šum okolo konštantného priemeru)

Vlasnosti bieleho šumu

  • šoky sú dočasné vo svojom efekte (mimo aktuálneho obdobia nemá žiaden vplyv), neovplyvní úroveň okolo priemeru, okolo ktorej sú generované budúce hodnoty

  • premenná vykazuje zmenu smeru k priemeru

  • šok nevplýva na prognózu

  • rozptyl chyby prognózy je ohraničený

Generovanie bieleho šumu

rnorm(100, 0, 1) = generovanie 100 náhodných hodnôt s priemerom 0 a odchýlkou 1

min(u) = kontrola pravidla 2σ, 3σ (tiež cez summary(u))

y <- 0.8 + u = generovanie bieleho šumu okolo konštantného priemeru

set.seed(12345)
u <- rnorm(100, 0, 1)

min(u)
[1] -2.380358
summary(u)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-2.3804 -0.5901  0.4837  0.2452  0.9004  2.4771 
plot(u, type = "l", col="coral"); abline(h = 0, lty = 2)

y <- 0.8 + u
plot(y, type = "l", col="seagreen"); abline(h = 0, lty = 2)

Funkcie reakcie na impulz a meranie časovej závislosti

Funkcia reakcie na impulz bieleho šumu je to, ako vplýva dnešný šok na premennú teraz, a ako na ňu bude vplývať v budúcnosti. TLDR; neexistuje vzťah medzi náhodnou zložkou ut a budúcimi hodnotami y v danom modeli. Šoky bieleho šumu sú v čase nekorelované, dnešný šok nemôže ovplyvniť budúce hodnoty y cez budúce hodnoty bieleho šumu. Proces bieleho šumu sa považuje za proces bez pamäti.

ACF = AutoCorrelation Function, zabudovaný z library(stats) sa nepoužíva, (vraj) je to kokotina

Z library(astsa)::acf2() máme dvojitý výstup - korelačná a parciálna korelačná funkcia (korelogram). Žiadna z hodnôt neprekročila hranicu, čo znamená, že hodnoty sú štatisticky nevýznamné

# acf(y, lag.max = 20)
# install.packages("astsa")
library(astsa)
acf2(y, max.lag=20) 

      [,1]  [,2]  [,3]  [,4] [,5]  [,6] [,7]  [,8]  [,9] [,10] [,11] [,12]
ACF  -0.02 -0.09 -0.10 -0.05 0.13 -0.02 0.08 -0.11 -0.11  0.01 -0.03  0.08
PACF -0.02 -0.09 -0.11 -0.07 0.11 -0.03 0.09 -0.09 -0.10 -0.02 -0.06  0.03
     [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20]
ACF   0.13 -0.19 -0.04 -0.13 -0.01  0.12  0.07 -0.01
PACF  0.15 -0.18 -0.01 -0.14 -0.09  0.06  0.08 -0.04

Ljung-Box test autokorelácie bieleho šumu okolo konštantného priemeru. (predpokladám, že AR teda nie je). Kedysi sa LB test dal robiť cez library(FitAR), ale už ho dali preč.

Box.test(y, type="Ljung-Box", lag=1)

    Box-Ljung test

data:  y
X-squared = 0.062381, df = 1, p-value = 0.8028
Box.test(y, type="Ljung-Box", lag=2)

    Box-Ljung test

data:  y
X-squared = 0.92803, df = 2, p-value = 0.6288

Kód, ktorý muž Lukáčikovej pacol do Teamsu, nastavenie Ljung-Box testu s Q-štatistikou. Nebude to fungovať, pretože LungBoxTest je súčasťou library(FitAR), ktorá už nie je dostupná. Ale pohoda. Máme skúsiť spojiť acf2(y, plot = F) a LjungBoxTest.

# LBQtab <- function(y, k=0, lag.max=16){
#  round(cbind("lag"=LjungBoxTest(y,k=k,lag.max=lag.max)[,1],
#              acf2(y,max.lag = lag.max,plot=F),
#              "Q-stat"=LjungBoxTest(y,k=k,lag.max=lag.max)[,2],
#              "p-value"=LjungBoxTest(y,k=k,lag.max=lag.max)[,3]),4)}
# LBQtab(y, lag.max=12, k=0)

Stacionarita a nestacionarita

Stacionarita

Pravdepodobnostné rozdelenie série náhodných premenných usporiadaných v čase, ktorá tvorí stochastický proces, zostáva rovnaké a nemení sa v čase. Časová závislosť medzi dvomi hodnotami nezávisí od období, ku ktorým prislúchajú, ale od ich vzdialenosti (nemenia sa s časom = sú časovo invariantné)

Ľudsky povedané, správanie premennej sa v čase nemení, len náhodne kolíše okolo nejakej úrovne. Teda:

  • priemer aj rozptyl sú konštantné

  • vzťahy medzi hodnotami (autokorelácia, autokovariancia) sa nemenia

Typy stacionarity

  • slabá (kovariančná) stacionarita: priemer, rozptyl a autokovarianca procesu musia byť časovo invariantné (nemenia sa s časom), je reštriktívnejšia

  • silná stacionarita: priemer, rozptyl, autokovarianca procesu A všetky momenty pravdepodobostného rozdelenia musia byť časovo invariantné (nemenia sa s časom)

Proces bieleho šumu je kovariančne stacionárny proces, pretože žiadna z jeho charakteristík nezávisí od času: E(yt) = β0, 𝛾0 = σ2, 𝛾j = 0 pre j > 0

Nestacionarita

Proces je nestacionárny, ak jedna z charakteristík procesu (priemer, rozptyl, autokovariancia) závisí od času (nie je časovo invariantná). Napríklad, rastúci trend znamená, že priemer nie je konštantný (teda priemer nie je časovo invariantný).

Príklady nestacionárnych procesov. Vo všetkých prípadoch sa s plynúcim časom premenná posúva ďalej od svojej počiatočnej hodnoty.

  • trendovo stacionárny proces: nestacionárny vzhľadom na svoj priemer

  • proces náhodnej prechádzky (random walk): nestacionárny vzhľadom na svoj rozptyl

  • proces náhodnej prechádzky s posunom: nestacionárny vzhľadom na svoj priemer aj rozptyl

Trendovo stacionárne procesy

Nestacionarita vzhľadom na priemer, yt fluktuuje okolo časového trendu t, nie okolo konštanty (priemeru).

Priemer sa mení s časom, nie je časovo invariantný, yt je nestacionárny proces. Krátkodobé fluktuácie bieleho šumu sú okolo trendu, detrendizáciou dostaneme z pôvodného procesu stacionárny rad (neanalyzujeme pôvodný proces, ale jeho odchýlky od trendu)

Postup detrendizácie: odhad modelu trendu, kde reziduály odhadu sú fluktuácie okolo trendu.

Platia vlastnosti procesu bieleho šumu okolo deterministického trendu.

  • šoky sú dočasné - mimo svojho obdobia nemá šok žiaden vplyv

  • priemer (deterministický trend) je hodnota, od ktorej sa premenná príliš nevzdiali

  • šok nevplýva na prognózu (najlepšou prognózou na ďalšie obdobie by mala byť ďalšia hodnota trendu)

  • rozptyl chyby prognózy (rozdiel medzi skutočnou a prognózovanou hodnotou) je ohraničnený (pravidlo 3σ)

  • autokorelačná funkcia by mala byť 0 (mimo aktuálneho obdobia) - ?

Proces náhodnej prechádzky

Nestacionarita vzhľadom na rozptyl

Proces náhodnej prechádzky s posunom (driftom)

tiridi bambambam

Ďalšie procesy vychádzajúce z bieleho šumu

Uvedené procesy sa využívajú pri prognózovaní, postup ich identifikácie a odhadu navrhli Box a Jenkins (podľa nich sa nazýva metodológia ARIMA modelov)

Procesy kĺzavých priemerov (moving averages)

Odlišujú sa podľa rádu oneskorenie bieleho šumu

Označenie procesu MA prvého rádu: MA(1) / ARMA(0, 1)

zapisuje sa s thetou (rovnička)

Stacionárne procesy, s krátkou (ohraničenou pamäťou), ich vlastnosti sa podobaju vlastnostiam bieleho šumu. Zvykne sa sledovať ich invertovateľnosť (ot autora - je to to iste abo suvisi s pojmom casovo invariantny?)

Autoregresné procesy (AR)

Odlišujú sa podľa rádu oneskorenia závislej premennej yt. Stacionárne procesy (podmienky)

Označenie autoregresného procesu prvého rádu: AR(1) / ARMA(1, 0)

Zapisuje sa s phí (rovnička)

Autoregresné procesy kĺzavých priemerov (ARMA)

Kombinácia dvoch predchádzajúcich typov procesov (AutoRegressive Moving Averages - ARMA). Odlišujú sa podľa:

  • rádu oneskorenia zaávislej premennej yt

  • rádu oneskorenia bieleho šumu

Označenie procesu prvého rádu: ARMA(1, 1)

Overuje sa ich stacionarita, invertovateľnosť (identifikácia pomocou ACF a PACF)

Integrované časové rady

Časový rad, ktorý nespĺňa podmienky stacionarity, ale pomoucou k-násobného diferencovania sa dá previesť na stacionárny rad

Označenie: yt ~ I(k)

Diferencovanie integrovaného radu yt prvého rádu = jeho stacionarizácia

Časové rady

Krátkodobé fluktuácie, abstrahovať od trendu a skúmať fluktuácie okolo pevne daného priemeru.

Integrované časové rady (?)

Príklad 1: unemployment rate

Časové rady unemp (miera nezamestnanosti) a GNP per capita z library(urca)::data(nporg). Z radov odstránime hodnoty NA.

Určenie typu procesu (stacionarita/nestacionarita). Asi začneme Ljung-Box testom pomocou jeho kódu LBQtab, ktorý nefunguje.

Následne autokorelačná (ACF) a parciálna autokorelačná funkcia (PACF) pomocou acf2(). Na základe tohto testu je to stacionárny rad

# LBQtab(unemp, lag.max=20, k=0)
acf2(unemp, max.lag=20) # na zaklade acf je to stacionarny rad, mohli by sme zacat analyzu Boxa a Jenkinsa modelom: 

     [,1]  [,2] [,3]  [,4] [,5] [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12] [,13]
ACF  0.86  0.65 0.48  0.31 0.18 0.09  0.01 -0.11 -0.19 -0.21 -0.22 -0.25 -0.24
PACF 0.86 -0.31 0.07 -0.18 0.05 0.03 -0.16 -0.18  0.08  0.04 -0.06 -0.16  0.09
     [,14] [,15] [,16] [,17] [,18] [,19] [,20]
ACF  -0.22 -0.22 -0.16  -0.1 -0.11 -0.12 -0.14
PACF -0.01 -0.03  0.10  -0.1 -0.16  0.03 -0.11
# ARIMA(2,0,4) - nieco s klzavymi priemermi - keby sme chceli (alebo SARIMA - package astsa) - chceme najlepis model na prognozovanie nezamestnanosti
sarima(unemp, 2, 0, 4) #, P, B, Q nemusime davat lebo su prednastavene (v helpe to vidis), 
initial  value 1.718423 
iter   2 value 1.343204
iter   3 value 1.262268
iter   4 value 1.130233
iter   5 value 0.999787
iter   6 value 0.976986
iter   7 value 0.969930
iter   8 value 0.964210
iter   9 value 0.963400
iter  10 value 0.963066
iter  11 value 0.962616
iter  12 value 0.962472
iter  13 value 0.962332
iter  14 value 0.962073
iter  15 value 0.961797
iter  16 value 0.961657
iter  17 value 0.961640
iter  18 value 0.961625
iter  19 value 0.961622
iter  20 value 0.961608
iter  21 value 0.961604
iter  22 value 0.961603
iter  23 value 0.961603
iter  24 value 0.961603
iter  25 value 0.961603
iter  26 value 0.961602
iter  27 value 0.961602
iter  28 value 0.961602
iter  28 value 0.961602
iter  28 value 0.961602
final  value 0.961602 
converged
initial  value 0.961391 
iter   2 value 0.961362
iter   3 value 0.960533
iter   4 value 0.959871
iter   5 value 0.959139
iter   6 value 0.955421
iter   7 value 0.952797
iter   8 value 0.952206
iter   9 value 0.949681
iter  10 value 0.946115
iter  11 value 0.945013
iter  12 value 0.943353
iter  13 value 0.943016
iter  14 value 0.942852
iter  15 value 0.942820
iter  16 value 0.942812
iter  17 value 0.942807
iter  18 value 0.942800
iter  19 value 0.942785
iter  20 value 0.942753
iter  21 value 0.942751
iter  22 value 0.942743
iter  23 value 0.942741
iter  24 value 0.942733
iter  25 value 0.942722
iter  26 value 0.942688
iter  27 value 0.942647
iter  28 value 0.942629
iter  29 value 0.942602
iter  30 value 0.942521
iter  31 value 0.942437
iter  32 value 0.942201
iter  33 value 0.942055
iter  34 value 0.941966
iter  35 value 0.941901
iter  36 value 0.941896
iter  37 value 0.941888
iter  38 value 0.941881
iter  39 value 0.941875
iter  40 value 0.941874
iter  41 value 0.941874
iter  42 value 0.941873
iter  43 value 0.941872
iter  44 value 0.941872
iter  45 value 0.941872
iter  45 value 0.941872
iter  45 value 0.941872
final  value 0.941872 
converged
<><><><><><><><><><><><><><>
 
Coefficients: 
      Estimate     SE t.value p.value
ar1     0.4799 0.2281  2.1039  0.0388
ar2    -0.4110 0.2122 -1.9369  0.0566
ma1     0.7206 0.2274  3.1689  0.0022
ma2     0.7968 0.1290  6.1776  0.0000
ma3     0.8811 0.1168  7.5456  0.0000
ma4     0.5504 0.1714  3.2117  0.0020
xmean   7.0301 1.1619  6.0507  0.0000

sigma^2 estimated as 6.259624 on 74 degrees of freedom 
 
AIC = 4.919153  AICc = 4.938094  BIC = 5.155642 
 

# sezonnost nezadavame lebo unemp nie je sezonny rad, vhodnejsie je rad unemp previest na time series (trz je to len vektor)
unemp <- ts(unemp, start=1890, frequency=1)
sarima(unemp, 2, 0, 4)
initial  value 1.718423 
iter   2 value 1.343204
iter   3 value 1.262268
iter   4 value 1.130233
iter   5 value 0.999787
iter   6 value 0.976986
iter   7 value 0.969930
iter   8 value 0.964210
iter   9 value 0.963400
iter  10 value 0.963066
iter  11 value 0.962616
iter  12 value 0.962472
iter  13 value 0.962332
iter  14 value 0.962073
iter  15 value 0.961797
iter  16 value 0.961657
iter  17 value 0.961640
iter  18 value 0.961625
iter  19 value 0.961622
iter  20 value 0.961608
iter  21 value 0.961604
iter  22 value 0.961603
iter  23 value 0.961603
iter  24 value 0.961603
iter  25 value 0.961603
iter  26 value 0.961602
iter  27 value 0.961602
iter  28 value 0.961602
iter  28 value 0.961602
iter  28 value 0.961602
final  value 0.961602 
converged
initial  value 0.961391 
iter   2 value 0.961362
iter   3 value 0.960533
iter   4 value 0.959871
iter   5 value 0.959139
iter   6 value 0.955421
iter   7 value 0.952797
iter   8 value 0.952206
iter   9 value 0.949681
iter  10 value 0.946115
iter  11 value 0.945013
iter  12 value 0.943353
iter  13 value 0.943016
iter  14 value 0.942852
iter  15 value 0.942820
iter  16 value 0.942812
iter  17 value 0.942807
iter  18 value 0.942800
iter  19 value 0.942785
iter  20 value 0.942753
iter  21 value 0.942751
iter  22 value 0.942743
iter  23 value 0.942741
iter  24 value 0.942733
iter  25 value 0.942722
iter  26 value 0.942688
iter  27 value 0.942647
iter  28 value 0.942629
iter  29 value 0.942602
iter  30 value 0.942521
iter  31 value 0.942437
iter  32 value 0.942201
iter  33 value 0.942055
iter  34 value 0.941966
iter  35 value 0.941901
iter  36 value 0.941896
iter  37 value 0.941888
iter  38 value 0.941881
iter  39 value 0.941875
iter  40 value 0.941874
iter  41 value 0.941874
iter  42 value 0.941873
iter  43 value 0.941872
iter  44 value 0.941872
iter  45 value 0.941872
iter  45 value 0.941872
iter  45 value 0.941872
final  value 0.941872 
converged
<><><><><><><><><><><><><><>
 
Coefficients: 
      Estimate     SE t.value p.value
ar1     0.4799 0.2281  2.1039  0.0388
ar2    -0.4110 0.2122 -1.9369  0.0566
ma1     0.7206 0.2274  3.1689  0.0022
ma2     0.7968 0.1290  6.1776  0.0000
ma3     0.8811 0.1168  7.5456  0.0000
ma4     0.5504 0.1714  3.2117  0.0020
xmean   7.0301 1.1619  6.0507  0.0000

sigma^2 estimated as 6.259624 on 74 degrees of freedom 
 
AIC = 4.919153  AICc = 4.938094  BIC = 5.155642

Príklad 2: IBM

Druhým pozorovaným časovým radom je rad akcií IBM.