Panelová data

Prurezová data:

\(y_i=\beta_0+\beta_1x_i+\epsilon_i\)

kde \(i\) predstavuje urcitou jednotku (firmu, stát, cloveka atd.), kterých máme N.

Casová rada

\(y_t=\beta_0+\beta_1x_t+\epsilon_t\)

Sledujeme konkrétní jednotku (firmu, stát, cloveka atd.) v case. Délka casové rady je T.

Panelová data

\(y_{it}=\beta_0+\beta_1x_{it}+\epsilon_{it}\)

Sledujeme jednotky \(i\) v case \(t\).

Nyní si predstavme, ze kazdá jednotka \(i\) obsahuje nejaký specifický faktor, oznacme jej jako \(c_i\). Tento faktor se nemeni v case, jinak by bylo \(c_{it}\).

\(y_{it}=\beta_0+\beta_1x_{it}+c_i+\epsilon_{it}\)

\(c_i\) predstavuje v case nemennou heterogenitu mezi promennými (individualni efekt). Napríklad pokud budeme zkoumat kriminalitu v rámci ruzných mest, tak \(c_i\) predstavuje nepozorovatelné vlastnosti daného mesta \(i\) a \(v_{it}\) zachycuje vsechny ostatni vlivy ovlivnujici \(y_{it}\), ktere jiz nejsme schopni vysvetlit.

Otázkou je, jak uchopit tento individuální efekt. V prípade, ze jej nezahrneme do odhadu, prejde \(c_i\) do nové náhodné slozky \(\epsilon_{it}\):

\(v_{it}=c_i+\epsilon_{it}\)

\(y_{it}=\beta_0+\beta_1x_{it}+v_{it}\)

Problémem nastane, kdyz bude existovat mezi \(c_i\) a \(x_{it}\) korelace \(Cor(c_i,x_{it})\neq 0\) Zároven platí, ze \(Cor(v_{it},x_{it})\neq 0\). Tím je porusen Gauss-Markov predpoklad o exogenite nezávislé promenné a odhad pomocí OLS bude nekonzistentní.

V takovém prípade musíme pouzít jinou metodu odhadu. Jako vhodné varianty muzeme pouzít:

V této cásti cvicení se budeme zaobírat fixed effect model.

Dalsí výhody panelových dat:

Rozlisujeme krátké a dlouhé panely

Síla testu je císlo mezi 0 a 1, které udává pravdepodobnost, ze pri neplatnosti nulové hypotézy dojde k jejímu zamítnutí, tedy pravdepodobnost odhalení neplatnosti nulové hypotézy. Platí, ze cím vyssí je síla testu, tím lépe.

install.packages("foreign")     # balik pro nahravani ruzných typu souboru STATA atd. 
install.packages("plm")         # balik pro panelová data
install.packages("lmtest")      # balík pro testy 
library(foreign)
library(plm)
library(lmtest)
library(sandwich)
getwd()
## [1] "C:/Users/Lukas/Desktop/Ekonometrie/AKM/NEW!!!/Cviceni/panel"
kravy=read.csv("C:/Users/Lukas/Desktop/Ekonometrie/AKM/NEW!!!/Cviceni/panel/kravy.csv")

Definujeme data jako panely.

kravy_panel=plm.data(kravy, index = c("FARM", "YEAR"))
attach(kravy_panel)

Balanced Panel 247 farms, 6 years (1993-1998), 1482 Observations Source: Author, Antonio Alvarez, Luis Arias

\(log(milk)=\beta_0+\beta_1 cows + \beta_2 land+ \beta_4 labor+\beta_5 feed+ \epsilon\)

OLS

ols=lm(log(MILK)~COWS+LAND+LABOR+FEED, data = kravy)
summary(ols)
## 
## Call:
## lm(formula = log(MILK) ~ COWS + LAND + LABOR + FEED, data = kravy)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.2781 -0.1390  0.0506  0.1880  0.7251 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.042e+01  2.480e-02 420.037  < 2e-16 ***
## COWS         3.751e-02  1.760e-03  21.306  < 2e-16 ***
## LAND        -3.100e-03  1.630e-03  -1.902   0.0574 .  
## LABOR        1.163e-01  1.524e-02   7.629 4.22e-14 ***
## FEED         3.084e-06  3.422e-07   9.011  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2676 on 1477 degrees of freedom
## Multiple R-squared:  0.8275, Adjusted R-squared:  0.827 
## F-statistic:  1771 on 4 and 1477 DF,  p-value: < 2.2e-16

Fixed effects model

\(y_{it}=\alpha+\beta_1x_{it}+c_i+v_{it}\)

V prípade existence endogenity máme moznost vyuzít jeden ze 3 mozných odhadu

  • LSDV regresi
  • first differences (FD)
  • fixed effects estimátor nebo také within estimátor (FE)

Fixed effects odhadnutý pomocí Least squares dummy variable model

V tomto prípade je kazdé jednotce \(i\) prirazena dummy promenná, ktera vlastne priradíme kazdému mestu dummy promennou.

\(y_{it}=\beta_1x_{it}+\alpha_i D_i+v_{it}\)

lsdv=lm(log(MILK)~COWS+LAND+LABOR+FEED+factor(FARM)-1, data=kravy_panel)
Box.test(lsdv$residua,lag = 10)   # test autokorelace
## 
##  Box-Pierce test
## 
## data:  lsdv$residua
## X-squared = 390.35, df = 10, p-value < 2.2e-16
bptest(lsdv, data = kravy_panel, studentize=F)
## 
##  Breusch-Pagan test
## 
## data:  lsdv
## BP = 1166.4, df = 250, p-value < 2.2e-16
robustni_odhad=coeftest(lsdv, vcov = vcovHAC)

Nevýhody tohoto odhadu nastávají pro prípad vyssích hodnot N. Ztrácíme stupne volnosti. \(N+k\) parametru.

Není urovnova konstanta viz. dummy trap. Pokud jsou v v modelu zastoupené promenné, které se nemení v case (rasa, pohlaví), tak musí být z odhadu vyrazeny z duvodu dokonalé multikolinearity.

Fixed effects - “Witin Estimator” nebo také Fixed effect estimator

Dalsí mozností je vyuzití fixed effect estimátoru, nebo téz Within estimatoru. Název je mozná lehce matoucí, ale FD, LSDV i FE estimátor patrí do skupiny fixed effect model.

Tato metoda je zalozena na tzv. “time-demeaned” (casove centrované) úprave, nebo také within transformation:

\(y_{it}-\bar{y_i}=\beta_1(x_{it}-\bar{x_i})+c_i-\bar{c_i}+v_{it}-\bar{v_i}\)

kdy

\(\bar{y_i}=\dfrac{\sum_{t=1}^{T}y_{it}}{T}\)

Jedna se o aritmeticky prumer pro prurezovou jednotky \(i\) pres vsechna casová období \(t\).

\(c_i-\bar{c_i}=0\)

Získáme tak:

\(\ddot{y}_{it}=\beta_1 \ddot{x}_{it}+\ddot\epsilon_{it}\)

a muzeme aplikovat OLS.

fix_wit = plm(log(MILK)~COWS+LAND+LABOR+FEED, data=kravy_panel, model="within")
summary(fix_wit)
## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = log(MILK) ~ COWS + LAND + LABOR + FEED, data = kravy_panel, 
##     model = "within")
## 
## Balanced Panel: n = 247, T = 6, N = 1482
## 
## Residuals:
##     Min.  1st Qu.   Median  3rd Qu.     Max. 
## -0.55500 -0.05540  0.00578  0.06190  0.47400 
## 
## Coefficients:
##          Estimate  Std. Error t-value Pr(>|t|)    
## COWS   3.2686e-02  1.6524e-03 19.7808   <2e-16 ***
## LAND   1.7403e-03  1.5817e-03  1.1002   0.2714    
## LABOR -1.4806e-02  2.1013e-02 -0.7046   0.4812    
## FEED   2.8306e-06  2.5307e-07 11.1853   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    49.745
## Residual Sum of Squares: 16.789
## R-Squared:      0.66251
## Adj. R-Squared: 0.59396
## F-statistic: 604.116 on 4 and 1231 DF, p-value: < 2.22e-16

\(\bf Breusch-Godfrey\) \(\bf test\) \(\bf pro\) \(\bf panelová\) \(\bf data\)

  • \(H_0\) znacní, ze zde není sériová korelace
pbgtest(fix_wit, order=1)
## 
##  Breusch-Godfrey/Wooldridge test for serial correlation in panel
##  models
## 
## data:  log(MILK) ~ COWS + LAND + LABOR + FEED
## chisq = 31.436, df = 1, p-value = 2.061e-08
## alternative hypothesis: serial correlation in idiosyncratic errors

Testování heteroskedasticity

\(\bf Breusch-Pagan~test\)

  • \(H_0\) homoskedasticita
bptest(fix_wit, data = kravy_panel, studentize=F)
## 
##  Breusch-Pagan test
## 
## data:  fix_wit
## BP = 200.59, df = 4, p-value < 2.2e-16

Nyní kdyz porovnáme výsledky z FE s výsledky z LSDV, muzeme videt, ze odhady parametru jsou shodné. Nevýhodou je, ze opet prijdeme o vliv promenných, které se nemení v case (pohlavi, rasa atd.)

Fixed effects - “first difference” FD

Tato metoda pracuje na principu odstranení nepozorované promenné \(a_i\) pomocí prvních diferencí.

\(\Delta y_{it}=\beta_1(x_{it}-x_{it-1})+\epsilon_{it}-\epsilon_{it-1}\)

Z rovnice muzeme videt, ze v prípade, kdy je jako regresor promenná, která se nemení v case, jako je napríklad pohlaví, nebo vzdelání, tak jsou tyto promenné vylouceny z regrese stejne jako \(a_i\). V prípade, kdy nás nezajímá daná promenná, to není az tak velký problém.

Opet pouzijeme funkci \(\textbf{plm}\) pouze zmeníme model z pooling na fd.

fix_fd = plm(log(MILK)~COWS+LAND+LABOR+FEED, data=kravy_panel, model="fd")
summary(fix_fd)
## Oneway (individual) effect First-Difference Model
## 
## Call:
## plm(formula = log(MILK) ~ COWS + LAND + LABOR + FEED, data = kravy_panel, 
##     model = "fd")
## 
## Balanced Panel: n = 247, T = 6, N = 1482
## Observations used in estimation: 1235
## 
## Residuals:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -0.5690 -0.0305  0.0293  0.0315  0.0975  0.6160 
## 
## Coefficients:
##         Estimate Std. Error t-value  Pr(>|t|)    
## COWS  2.8537e-02 1.5064e-03 18.9440 < 2.2e-16 ***
## LAND  3.9670e-03 1.3765e-03  2.8819  0.004021 ** 
## LABOR 9.8000e-03 1.8812e-02  0.5209  0.602510    
## FEED  2.3768e-06 2.3685e-07 10.0352 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    23.57
## Residual Sum of Squares: 16.276
## R-Squared:      0.37182
## Adj. R-Squared: 0.37029
## F-statistic: 183.899 on 3 and 1231 DF, p-value: < 2.22e-16
pbgtest(fix_fd, order=1)
## 
##  Breusch-Godfrey/Wooldridge test for serial correlation in panel
##  models
## 
## data:  log(MILK) ~ COWS + LAND + LABOR + FEED
## chisq = 3.2595, df = 1, p-value = 0.07101
## alternative hypothesis: serial correlation in idiosyncratic errors
bptest(fix_fd, data = kravy_panel, studentize=F)
## 
##  Breusch-Pagan test
## 
## data:  fix_fd
## BP = 200.59, df = 4, p-value < 2.2e-16

Random effects

\(y_{it}=\alpha+x_{it}\beta_1+c_i+\epsilon_{it}\)

kdy \(c_i\) predstavuje random heterogenity pro jednotku \(i\). Získáme náhodnou slozku \(v_{it}=c_i+\epsilon_{it}\)

Zásadní predpoklad je, ze \(Cov(x_{it},v_{it})=0\). Tedy opak toho co predpokládáme u fix effect model.

\(\bf Výhody\)

  • muzeme zkoumat v case nemenné promenní (pohlaví, rasa)
  • promenné s malou variabilitou
  • snízíme mnozství odhadovaných parametru

\(\bf Nevýhody\)

  • moznost nekonzistetního odhadu, pokud není splneno:

\(Cov(x_{it},v_{it})=0\).

Pokud je vsak tento predpoklad splnen, muzeme pouzít OLS, která poskytne nezkreslené a konzistentní odhady. Ne vsak vydatné! Problémem je prítomnost autokorelace v \(v_{it}\)

Pro odhad je nutné pouzít (F)GLS. Jedná se o tzv. casového kvazicentrování (quasi-time demeaning).

\(cor(v_{it},v_{is})=(\dfrac{\sigma^2_u}{\sigma^2_u+\sigma^2_{v}})\)

\(y_{it}-\lambda \bar{y}_i =(y_{it}-\lambda \bar{y}_i)\beta+(v_{it}-\lambda \bar{v_i})\)

kdy

\(\lambda=\dfrac{1}{1+T(\sigma^2_u/\sigma^2_{v})^(0.5)}\)

A toto je \(\textbf{RANDOM EFFECT ESTIMATOR}\)

Ve vetsine statistických programech je vsak tento odhad jiz naimplementován.

rem =plm(log(MILK)~COWS+LAND+LABOR+FEED, data=kravy_panel, model="random")
summary(rem)
## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = log(MILK) ~ COWS + LAND + LABOR + FEED, data = kravy_panel, 
##     model = "random")
## 
## Balanced Panel: n = 247, T = 6, N = 1482
## 
## Effects:
##                   var std.dev share
## idiosyncratic 0.01364 0.11678  0.19
## individual    0.05820 0.24124  0.81
## theta: 0.8061
## 
## Residuals:
##    Min. 1st Qu.  Median 3rd Qu.    Max. 
## -0.6270 -0.0617  0.0121  0.0761  0.4720 
## 
## Coefficients:
##               Estimate Std. Error  t-value Pr(>|t|)    
## (Intercept) 1.0560e+01 3.6616e-02 288.3986  < 2e-16 ***
## COWS        3.5215e-02 1.4874e-03  23.6760  < 2e-16 ***
## LAND        1.1965e-03 1.4722e-03   0.8128  0.41648    
## LABOR       3.8866e-02 1.7845e-02   2.1780  0.02957 *  
## FEED        2.7300e-06 2.4313e-07  11.2286  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    70.923
## Residual Sum of Squares: 20.508
## R-Squared:      0.71085
## Adj. R-Squared: 0.71007
## F-statistic: 907.765 on 4 and 1477 DF, p-value: < 2.22e-16

Který model pouzít?

Jak se máme rozhodnout, který model pouzít?

  • Pooled vs. Random \(\rightarrow\) LM test
  • Fixed vs. Pool \(\rightarrow\) Chow test
  • Fixed vs. Random \(\rightarrow\) Hausman test

Fixed model vs. Random model

Hausman test

V prípade, ze \(Cor(v_{it},x_{it})= 0\) budou odhady pomocí FE, respektive RE velmi podobné. Pokud vsak bude platit, ze \(Cor(v_{it},x_{it})\neq 0\) budou rozdíly mezi odhady statisticky významé. V tomto prípade také zálezí na N a T.

Metoda \(H_0:Exogenní\) \(H_1:Endogenní\)
FE Konzistetní a nevydatný Konzistetní
RE Konzistetní a vydatný Nekonzistetní
phtest(fix_wit, rem)   # predpoklad homoskedasticity a NEautokorelace
## 
##  Hausman Test
## 
## data:  log(MILK) ~ COWS + LAND + LABOR + FEED
## chisq = 32.076, df = 4, p-value = 1.846e-06
## alternative hypothesis: one model is inconsistent

Ruzné robustní verze Hausmanova testu

phtest(form, data = kravy_panel, method = "aux",vcov = function(x) vcovHC(x, method="arellano", type="HC3"))

Podle eko. teorie

\(log(milk)=\beta_0+\beta_1 log(cows) + \beta_2 log(land) + \beta_4 log(labor)+\beta_5 log(feed)+ \epsilon\)

Lfix_wit = plm(log(MILK)~log(COWS)+log(LAND)+log(LABOR)+log(FEED), data=kravy_panel, model="within")
summary(Lfix_wit)
## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = log(MILK) ~ log(COWS) + log(LAND) + log(LABOR) + 
##     log(FEED), data = kravy_panel, model = "within")
## 
## Balanced Panel: n = 247, T = 6, N = 1482
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -0.467000 -0.041300  0.000129  0.045200  0.370000 
## 
## Coefficients:
##            Estimate Std. Error t-value Pr(>|t|)    
## log(COWS)  0.662001   0.024678 26.8251  < 2e-16 ***
## log(LAND)  0.037352   0.016133  2.3153  0.02076 *  
## log(LABOR) 0.030400   0.023208  1.3099  0.19048    
## log(FEED)  0.382510   0.012017 31.8310  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    49.745
## Residual Sum of Squares: 8.1611
## R-Squared:      0.83594
## Adj. R-Squared: 0.80262
## F-statistic: 1568.11 on 4 and 1231 DF, p-value: < 2.22e-16

Nonparametric robust covariance matrix estimators a la Driscoll and Kraay for panel models with cross-sectional and serial correlation.

coeftest(Lfix_wit, vcov=vcovSCC)
## 
## t test of coefficients:
## 
##             Estimate Std. Error t value  Pr(>|t|)    
## log(COWS)  0.6620010  0.0070339 94.1163 < 2.2e-16 ***
## log(LAND)  0.0373525  0.0071010  5.2602 1.696e-07 ***
## log(LABOR) 0.0303995  0.0114843  2.6470  0.008224 ** 
## log(FEED)  0.3825104  0.0256920 14.8883 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(Lfix_wit, vcovSCC(Lfix_wit, type="HC1", maxlag=4))
## 
## t test of coefficients:
## 
##             Estimate Std. Error  t value  Pr(>|t|)    
## log(COWS)  0.6620010  0.0048625 136.1441 < 2.2e-16 ***
## log(LAND)  0.0373525  0.0045501   8.2092 5.566e-16 ***
## log(LABOR) 0.0303995  0.0079653   3.8165 0.0001421 ***
## log(FEED)  0.3825104  0.0165189  23.1559 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(Lfix_wit, vcovHC(Lfix_wit, method = "arellano")) 
## 
## t test of coefficients:
## 
##            Estimate Std. Error t value Pr(>|t|)    
## log(COWS)  0.662001   0.034062 19.4352  < 2e-16 ***
## log(LAND)  0.037352   0.017088  2.1859  0.02901 *  
## log(LABOR) 0.030400   0.024241  1.2541  0.21005    
## log(FEED)  0.382510   0.017235 22.1939  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Testování modelu na splnení GM

Testování sériové korelace

\(\bf Breusch-Godfrey\) \(\bf test\) \(\bf pro\) \(\bf panelová\) \(\bf data\)

  • \(H_0\) znacní, ze zde není sériová korelace
pbgtest(fix_wit, order=1)
## 
##  Breusch-Godfrey/Wooldridge test for serial correlation in panel
##  models
## 
## data:  log(MILK) ~ COWS + LAND + LABOR + FEED
## chisq = 31.436, df = 1, p-value = 2.061e-08
## alternative hypothesis: serial correlation in idiosyncratic errors
pbgtest(rem, order=1)
## 
##  Breusch-Godfrey/Wooldridge test for serial correlation in panel
##  models
## 
## data:  log(MILK) ~ COWS + LAND + LABOR + FEED
## chisq = 123.58, df = 1, p-value < 2.2e-16
## alternative hypothesis: serial correlation in idiosyncratic errors

Testování heteroskedasticity

\(\bf Breusch-Pagan test\)

  • \(H_0\) homoskedasticita
bptest(fix_wit, data = Crime_panel, studentize=F)
## 
##  Breusch-Pagan test
## 
## data:  fix_wit
## BP = 200.59, df = 4, p-value < 2.2e-16
bptest(crmrte ~ polpc++prbarr+prbconv+prbpris, data = Crime.set, studentize=F)

Testování cross-sectional dependency

Modely, které jsme si zde predstavily jsou konstruovány spíse pro mikroekonomickou analýzu, kdy predpokládáme, ze prvotní výber z populace je náhodný. Problém nastává v prípade makroekonomických dat. Predstavte si, ze chcete zkoumat napríklad

\begin{equation} HDP_{it}=\beta_0+\beta_1 K_{it}+\beta_2 E_{it}+\epsilon_{it} \end{equation}

Je zrejmé, ze prvnotní výber není náhodný. Uvedený vztah chceme analyzovat napríklad pro zeme EU. Z tohoto duvodu je velmi pravdepodobné, ze mezi zememi bude existovat vztah, který se projeví v náhodné slozce.

Problém existence prurezové závislosti v panelových datech je relativne nová oblast a patrí spíse do pokrocilých kurzu. Pro vás je vsak dulezité, ze v prípade existence prurezové závislosti mohou vést výse uvedené metody k \({\bf NEKONZISTETNÍMU}\) odhadum.

Pro testování prurezové závislosti muzeme vyuzít napríklad CD (Cross-dependence) test podle Pesaran.

  • \(H_0:\) není prurezová závislost
pcdtest(rem, test = c("cd"))

Testování stacionarity

Tato cást patrí pod casové rady, které jeste neznáte. Pro vás je to spíse z duvodu, na co si dát pozor. Casové rady mezi jsou relativne oblíbené (HDP, ceny aktiv, atd.). Bohuzel jednou z nejcastejsích chyb je aplikace regresní analýzy na tzv. nestacionární data. V takovém prípade jsou výsledky mylné.

library(tseries)

Dickey-Fuller test

  • \(H_0:\) casová rada není stacionární
adf.test(kravy_panel$MILK, k=2)

V prípade existence jednotkového korene, je nutné data tzv. stacionarizovat. Stacionarizaci provedeme pomocí první diference.

\begin{equation} \Delta y_{it}=y_{it}-y_{it-1} \end{equation}

Logaritmická diference

Co delat v prípade autokorelace a heteroskedasticity

Robust covariance matrix estimation (Sandwich estimator)

Funkce pro odhad robustní kovariancní matice je vcovHC Dále existuje více robustních odhadu, kdy kazdý má svá urcitá specifika a proto je vhodnejsí pro ten, ci onen prípad:

  • “white1” - je robustní pro heteroskedasticitu, ne vsak pro autokorelaci vhodný pro random effects.
  • “arellano” - robustní jak pro heteroskedasticitu, tak autokorelaci, vhodný pro fixed effects.

  • HC0 - heteroskedasticity consistent. The default.
  • HC1,HC2, HC3 - Recommended for small samples. HC3 gives less weight to influential observations.
  • HC4 -small samples with influential observations
  • HAC - heteroskedasticity and autocorrelation consistent

Heteroskedasticity: Fixed effects

Heteroskedasticity consistent coefficients

coeftest(fix_wit) 
coeftest(fix_wit, vcovHC) 
coeftest(fix_wit, vcovHC(fix_wit, type = "HC3")) 

Heteroskedasticity: Random effects

coeftest(rem) 
coeftest(rem, vcovHC) # Heteroskedasticity consistent standard errors
coeftest(rem, vcovHC(rem, method = "arellano", cluester = "FARM")) 
coeftest(rem, vcovHC(rem, type = "HC3")) 

Dodatecný príklad

wage=read.dta("C:/Users/Lukas/Desktop/Ekonometrie/AKM/cviceni 2016/cviceni 9/wagepan.dta")
library("plm")
wage_panel=plm.data(wage, index = c("nr", "year"))
## Warning: use of 'plm.data' is discouraged, better use 'pdata.frame' instead
regrese_pool =plm(lwage~educ+black+hisp+exper+expersq+married+union, data=wage_panel, model="pooling")
summary(regrese_pool)
## Pooling Model
## 
## Call:
## plm(formula = lwage ~ educ + black + hisp + exper + expersq + 
##     married + union, data = wage_panel, model = "pooling")
## 
## Balanced Panel: n = 545, T = 8, N = 4360
## 
## Residuals:
##    Min. 1st Qu.  Median 3rd Qu.    Max. 
## -5.2700 -0.2490  0.0332  0.2960  2.5600 
## 
## Coefficients:
##                Estimate  Std. Error t-value  Pr(>|t|)    
## (Intercept) -0.03470560  0.06456900 -0.5375    0.5910    
## educ         0.09938779  0.00467760 21.2476 < 2.2e-16 ***
## black       -0.14384172  0.02355950 -6.1055 1.114e-09 ***
## hisp         0.01569798  0.02081119  0.7543    0.4507    
## exper        0.08917906  0.01011105  8.8200 < 2.2e-16 ***
## expersq     -0.00284866  0.00070736 -4.0272 5.742e-05 ***
## married      0.10766559  0.01569647  6.8592 7.897e-12 ***
## union        0.18007255  0.01712053 10.5179 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    1236.5
## Residual Sum of Squares: 1005.8
## R-Squared:      0.18659
## Adj. R-Squared: 0.18528
## F-statistic: 142.613 on 7 and 4352 DF, p-value: < 2.22e-16
regrese_fix=plm(lwage~educ+black+hisp+exper+expersq+married+union, data=wage_panel, model="within")
summary(regrese_fix)
## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = lwage ~ educ + black + hisp + exper + expersq + 
##     married + union, data = wage_panel, model = "within")
## 
## Balanced Panel: n = 545, T = 8, N = 4360
## 
## Residuals:
##     Min.  1st Qu.   Median  3rd Qu.     Max. 
## -4.17000 -0.12600  0.00925  0.16000  1.47000 
## 
## Coefficients:
##            Estimate  Std. Error t-value  Pr(>|t|)    
## exper    0.11684669  0.00841968 13.8778 < 2.2e-16 ***
## expersq -0.00430089  0.00060527 -7.1057 1.422e-12 ***
## married  0.04530333  0.01830968  2.4743   0.01339 *  
## union    0.08208713  0.01929073  4.2553 2.138e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    572.05
## Residual Sum of Squares: 470.2
## R-Squared:      0.17804
## Adj. R-Squared: 0.059852
## F-statistic: 206.375 on 4 and 3811 DF, p-value: < 2.22e-16

Pokud chceme zjistit, jaký vliv má vzdelání na mzdu, tak zde máme smulu.

regrese_random=plm(lwage~educ+black+hisp+exper+expersq+married+union, data=wage_panel, model="random")
summary(regrese_random)
## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = lwage ~ educ + black + hisp + exper + expersq + 
##     married + union, data = wage_panel, model = "random")
## 
## Balanced Panel: n = 545, T = 8, N = 4360
## 
## Effects:
##                  var std.dev share
## idiosyncratic 0.1234  0.3513 0.539
## individual    0.1053  0.3246 0.461
## theta: 0.6426
## 
## Residuals:
##    Min. 1st Qu.  Median 3rd Qu.    Max. 
## -4.5800 -0.1450  0.0235  0.1860  1.5400 
## 
## Coefficients:
##                Estimate  Std. Error t-value  Pr(>|t|)    
## (Intercept) -0.10746420  0.11070573 -0.9707 0.3317420    
## educ         0.10122461  0.00891329 11.3566 < 2.2e-16 ***
## black       -0.14413069  0.04761483 -3.0270 0.0024843 ** 
## hisp         0.02015107  0.04260112  0.4730 0.6362245    
## exper        0.11211949  0.00826087 13.5724 < 2.2e-16 ***
## expersq     -0.00406885  0.00059183 -6.8751 7.075e-12 ***
## married      0.06279512  0.01677285  3.7439 0.0001836 ***
## union        0.10737885  0.01783001  6.0224 1.861e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    656.91
## Residual Sum of Squares: 539.82
## R-Squared:      0.17824
## Adj. R-Squared: 0.17692
## F-statistic: 134.85 on 7 and 4352 DF, p-value: < 2.22e-16
pbgtest(regrese_fix, order=1)    #test autokorelace 
## 
##  Breusch-Godfrey/Wooldridge test for serial correlation in panel
##  models
## 
## data:  lwage ~ educ + black + hisp + exper + expersq + married + union
## chisq = 6.9358, df = 1, p-value = 0.008449
## alternative hypothesis: serial correlation in idiosyncratic errors
pbgtest(regrese_random, order=1)
## 
##  Breusch-Godfrey/Wooldridge test for serial correlation in panel
##  models
## 
## data:  lwage ~ educ + black + hisp + exper + expersq + married + union
## chisq = 113.42, df = 1, p-value < 2.2e-16
## alternative hypothesis: serial correlation in idiosyncratic errors
form <- lwage~educ+black+hisp+exper+expersq+married+union
phtest(form, data = wage_panel, method = "aux", vcov = function(x) vcovHC(x, method="arellano")) # Robustní Hausman test
## 
##  Regression-based Hausman test, vcov: function(x) vcovHC(x, method
##  = "arellano")
## 
## data:  form
## chisq = 81.83, df = 4, p-value < 2.2e-16
## alternative hypothesis: one model is inconsistent

Priklad

getwd()
## [1] "C:/Users/Lukas/Desktop/Ekonometrie/AKM/NEW!!!/Cviceni/panel"
load("C:/Users/Lukas/Desktop/Ekonometrie/AKM/NEW!!!/Cviceni/panel/crime4.RData")

Zadefinujeme jako panel

crime_panel=plm.data(data, index = c("county", "year"))
attach(crime_panel)

\(log(crmrte)=\beta_0+\beta_1 log(prbarr)+\beta_2 log(prbconv)+ \beta_3 polpc+\beta_4 log(wmfg)+\epsilon\)

ols=lm(log(crmrte)~I(log(prbarr))+I(log(prbconv))+polpc+I(log(wmfg)))
summary(ols)
## 
## Call:
## lm(formula = log(crmrte) ~ I(log(prbarr)) + I(log(prbconv)) + 
##     polpc + I(log(wmfg)))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.66161 -0.21828  0.03972  0.24348  1.18122 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -6.49514    0.31264 -20.775  < 2e-16 ***
## I(log(prbarr))  -0.74782    0.03784 -19.765  < 2e-16 ***
## I(log(prbconv)) -0.61161    0.02862 -21.372  < 2e-16 ***
## polpc           71.51029    6.35513  11.252  < 2e-16 ***
## I(log(wmfg))     0.24436    0.05669   4.311 1.89e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3797 on 625 degrees of freedom
## Multiple R-squared:  0.5633, Adjusted R-squared:  0.5605 
## F-statistic: 201.5 on 4 and 625 DF,  p-value: < 2.2e-16

Jelikoz máme panelová data sestavíme fixed effect model

\(log(crmrte)_{it}=\beta_0+\beta_1 log(prbarr)_{it}+\beta_2 log(prbconv)_{it}+ \beta_3 polpc_{it}+\beta_4 log(wmfg)_{it}+c_i+\epsilon_{it}\)

kde \(u_i\) predstavuje nepozorovanou heterogenitu daného (mesta) county.

Klidne si muzeme predstavit, ze platí:

\(c_i=Z_i'\alpha+e_i\)

Kde \(Z\) predstavuje nepozorovane faktory pro dané mesto (prumysl, poloha, velikost mesta) a \(e\sim iid(0,\sigma ^2)\).

\(log(crmrte)_{it}=\beta_0+\beta_1 log(prbarr)_{it}+\beta_2 log(prbconv)_{it}+ \beta_3 polpc_{it}+\beta_4 log(wmfg)_{it}+v_{it}\)

Kdy \(v_{it}=c_i+\epsilon_{it}\) je composite error.

Muzeme predpokádat, ze nekterý z faktoru v \(Z\) bude korelován s poctem policistu, napr. cím vetsí mesto, tím je potreba více policistu.

\(Cor(Z_i,polpc_{it})\neq0\implies Cor(c_i,polpc_{it})\neq0 \implies Cor(v_{it},polpc_{it})\neq0\)

Je porusena podmínka exogenity \(E(v_{it}|X)\neq 0 ~~ \forall i,t\)

Odhad FE modelu

Jako první pouzíjeme FD estimátor

FD=plm(log(crmrte)~I(log(prbarr))+I(log(prbconv))+polpc+I(log(wmfg)),data = crime_panel, model = "fd")
summary(FD)
## Oneway (individual) effect First-Difference Model
## 
## Call:
## plm(formula = log(crmrte) ~ I(log(prbarr)) + I(log(prbconv)) + 
##     polpc + I(log(wmfg)), data = crime_panel, model = "fd")
## 
## Balanced Panel: n = 90, T = 7, N = 630
## Observations used in estimation: 540
## 
## Residuals:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -0.86600 -0.08280  0.00834  0.00307  0.08990  0.93800 
## 
## Coefficients:
##                  Estimate Std. Error  t-value Pr(>|t|)    
## I(log(prbarr))  -0.301331   0.031954  -9.4302   <2e-16 ***
## I(log(prbconv)) -0.212107   0.019411 -10.9274   <2e-16 ***
## polpc           49.760513   4.574906  10.8768   <2e-16 ***
## I(log(wmfg))    -0.037628   0.084575  -0.4449   0.6566    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    22.197
## Residual Sum of Squares: 16.081
## R-Squared:      0.27577
## Adj. R-Squared: 0.27171
## F-statistic: 67.9526 on 3 and 536 DF, p-value: < 2.22e-16

Dále pouzijeme FE estimátor

FE=plm(log(crmrte)~I(log(prbarr))+I(log(prbconv))+polpc+I(log(wmfg)),data = crime_panel, model = "within")
summary(FE)
## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = log(crmrte) ~ I(log(prbarr)) + I(log(prbconv)) + 
##     polpc + I(log(wmfg)), data = crime_panel, model = "within")
## 
## Balanced Panel: n = 90, T = 7, N = 630
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -0.813000 -0.077500  0.000141  0.075200  0.783000 
## 
## Coefficients:
##                  Estimate Std. Error  t-value  Pr(>|t|)    
## I(log(prbarr))  -0.336847   0.033306 -10.1137 < 2.2e-16 ***
## I(log(prbconv)) -0.265954   0.021589 -12.3189 < 2.2e-16 ***
## polpc           55.942371   4.088919  13.6815 < 2.2e-16 ***
## I(log(wmfg))    -0.131876   0.047928  -2.7515  0.006133 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    17.991
## Residual Sum of Squares: 12.272
## R-Squared:      0.31788
## Adj. R-Squared: 0.19953
## F-statistic: 62.4473 on 4 and 536 DF, p-value: < 2.22e-16

Odhad parametru \(\beta_3\) je stále kladný. O cem to svedcí? Nás predpoklad o modelu FE nebude správný. Prícinou muze být napr. prítomnost prurezové závislosti

TOTO JEZ NENÍ NÁPLNÍ KURZU!!!!

Ale ukázeme si, jak za urcitých okolností muze být nekonzistntní odhad FE zpusopben tím, ze výber není náhodný a mezi panely existuje prurezová závislost.

Otestujeme na purezovou závislosti

pcdtest(FE, test = c("cd"))
## 
##  Pesaran CD test for cross-sectional dependence in panels
## 
## data:  log(crmrte) ~ I(log(prbarr)) + I(log(prbconv)) + polpc + I(log(wmfg))
## z = 27.596, p-value < 2.2e-16
## alternative hypothesis: cross-sectional dependence

Zamítáme H0 ve prospech H1, která svedcí pro data zatízena prurezovou závislostí.

Proto provedeme odhad robustní vuci prurezové závislosti tzv. CCE estimátor viz Pesaran (2006)

cce=pcce(log(crmrte)~I(log(prbarr))+I(log(prbconv))+polpc+I(log(wmfg)),data = crime_panel, model = "mg")
summary(cce)
## Common Correlated Effects model
## Call:
## pcce(formula = log(crmrte) ~ I(log(prbarr)) + I(log(prbconv)) + 
##     polpc + I(log(wmfg)), data = crime_panel, model = "mg")
## 
## Balanced Panel: n = 90, T = 7, N = 630
## 
## Residuals:
##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## -7.188e-11 -9.579e-13  1.109e-14 -2.822e-14  9.618e-13  1.274e-10 
## 
## Coefficients:
##                    Estimate  Std. Error z-value  Pr(>|z|)    
## I(log(prbarr))  -0.20518613  0.05120257 -4.0073 6.141e-05 ***
## I(log(prbconv)) -0.06997564  0.04333262 -1.6148    0.1063    
## polpc            0.00018714  0.00024394  0.7672    0.4430    
## I(log(wmfg))    -0.02497080  0.01607788 -1.5531    0.1204    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Total Sum of Squares: 206.38
## Residual Sum of Squares: 6.2004e-20
## HPY R-squared: 1

Estimátor CCE je robustní jak vuci autokorelaci, tak heteroskedasticite, takze odhadnuté std. error muzeme brát jako relevantní. Vidíme, ze parametr \(\beta_3\) jiz není statisticky významný a jeho hodnota je témer nulová. Toto byl motivacní príklad, máme zde jeste problém krátných casových rad, netestoval jsem stacionaritu, kointegraci atd. Cílem vsak bylo ukázat, jak dulezité jsou predpoklady o modelu. Je zrejmé, ze predpoklady FE modelu byly v tomto prípade nedostatecné.