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:

najít vhodnou(é) instrumentální promenné
pro nepozorované \(c_i\) najít vhodnou proxy promennou
vyuzít simultánní rovnice
vyuzít fixed effect model

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

Dalsí výhody panelových dat:

zvýsíme pocet stupnu volnosti a tím i sílu testu
moznost komplexnejsích modelu
u makroekonomických modelu krátké casové rady, mozné resení panely (problém je, ze metody odhadu, které se zde ucíme jsou pro mikroekonomická data, vysoké N a nízké T.

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)

Priklad 1

Budeme testovat, zdali má nezamestnanost vliv na kriminalitu.

getwd()

## [1] "C:/Users/Lukas/Desktop/Ekonometrie/AKM/NEW!!!/Cviceni/panel"

crime=read.csv("C:/Users/Lukas/Desktop/Ekonometrie/AKM/NEW!!!/Cviceni/panel/crime.csv")

Nejprve odhadneme model jako, tedy klasicky bereme data jako prurezová.

\(crmrte=\beta_0+\beta_1 unemp+\epsilon\)

attach(crime)

## The following object is masked _by_ .GlobalEnv:
## 
##     crime

ols=lm(crmrte~unemp)
summary(ols)

## 
## Call:
## lm(formula = crmrte ~ unemp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -51.686 -23.889  -7.961  17.522  78.297 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 103.2434     8.0587   12.81   <2e-16 ***
## unemp        -0.3077     0.9317   -0.33    0.742    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 29.99 on 90 degrees of freedom
## Multiple R-squared:  0.00121,    Adjusted R-squared:  -0.009888 
## F-statistic: 0.109 on 1 and 90 DF,  p-value: 0.742

Pro jednoduchost predpokladejme, ze jsou splneny GM. Tedy vychází nám, ze nezamestnanost nemá vliv kriminalitu a odhad parametru je záporný.

Nyní vezmeme v potaz cas. Bereme tak, ze data pocházejí z pooled independent crosss-section (jedná se samozrejme o panel, ale chci ukázat tu výstavbu).

\(crmrte=\beta_0+\gamma d_{87}+\beta_1 unemp+\epsilon\)

ols_cas=lm(crmrte~unemp+factor(year))

bptest(ols_cas)

## 
##  studentized Breusch-Pagan test
## 
## data:  ols_cas
## BP = 6.6535, df = 2, p-value = 0.03591

coeftest(ols_cas,vcovHC)

## 
## t test of coefficients:
## 
##                Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)    93.42026   11.05189  8.4529 5.082e-13 ***
## unemp           0.42655    1.06696  0.3998    0.6903    
## factor(year)87  7.94041    7.26146  1.0935    0.2771    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Jelikoz zamítáme na hladine významnosti 5\(\%\) nulovou hypotezu o homoskedasticite, tak provedeme robustni odhad chyb. Co nám výjde. Nezamestnanost stále nemá stat. významný vliv na kriinalitu, ale jiz máme kladné znaménko.

Intercept nám ríká, ze v roce by byla prumerná kriminalita pri unemp=0 93,4 testných cinu na 100 000 obyvatel. Behem 5 let doslo ke zmenám a vlivem casu vzrostla kriminalita o 7,9 cinu. Parametr vsak není stat. významný.

Co kdyz máme nejaký nepozorovaný faktor \(c_i\), nazveme ho faktorem mesta, který bude korelován s nezamestnaností, napr. poloha atd.

\(crmrte=\beta_0+\gamma d_{87}+\beta_1 unemp+c_i+\epsilon\)

Jelikoz máme panelová data, budeme predpokládat FE model a zdurazníme, ze se jedná o panely

\(crmrte_{it}=\beta_0+\gamma d_{87}+\beta_1 unemp_{it}+c_i+\epsilon_{it}\)

Nejprve Rku rekneme, ze máme panel. Prostorová specifikace je state a casová year.

crime_panel=pdata.frame(crime,index=c("state","year"))

LSDV

Jako prvni estimator FE modelu vyuzijeme LSDV

\(crmrte_{it}=\beta_0+\gamma d_{87}+\beta_1 unemp_{it}+\alpha_i D_i+\epsilon_{it}\)

lsdv=lm(crmrte~unemp+factor(state)+factor(year))

bptest(lsdv)

## 
##  studentized Breusch-Pagan test
## 
## data:  lsdv
## BP = 92, df = 47, p-value = 9.589e-05

coeftest(lsdv,vcovHC)

## 
## t test of coefficients:
## 
##                   Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)      51.489251  12.787185  4.0266 0.0002201 ***
## unemp             2.217996   1.209859  1.8333 0.0735334 .  
## factor(state)2   17.291718  11.554649  1.4965 0.1416585    
## factor(state)3    4.688526  13.127918  0.3571 0.7226931    
## factor(state)4    7.006747   8.961452  0.7819 0.4384746    
## factor(state)5   24.497990  10.657633  2.2986 0.0263342 *  
## factor(state)6   43.565590  18.540159  2.3498 0.0233348 *  
## factor(state)7    2.546190   8.842502  0.2879 0.7747377    
## factor(state)8   12.027173  22.492610  0.5347 0.5955378    
## factor(state)9  -10.251101  26.628720 -0.3850 0.7021186    
## factor(state)10  14.602237  10.833945  1.3478 0.1846144    
## factor(state)11  -2.070378  11.520312 -0.1797 0.8582011    
## factor(state)12  28.453762  28.231933  1.0079 0.3190335    
## factor(state)13  94.955124   8.210159 11.5656 6.251e-15 ***
## factor(state)14   8.476689   8.339076  1.0165 0.3149474    
## factor(state)15  32.483724   8.960375  3.6253 0.0007453 ***
## factor(state)16  69.437887   9.576927  7.2505 4.916e-09 ***
## factor(state)17  72.723821  30.003493  2.4238 0.0195359 *  
## factor(state)18  18.495352  16.143738  1.1457 0.2581259    
## factor(state)19  77.458579  24.442651  3.1690 0.0027821 ** 
## factor(state)20  67.753279  15.995503  4.2358 0.0001144 ***
## factor(state)21  -8.239871  14.835726 -0.5554 0.5814291    
## factor(state)22  -6.837449  12.704989 -0.5382 0.5931715    
## factor(state)23   8.183490  11.847508  0.6907 0.4933601    
## factor(state)24  37.096057  14.218266  2.6090 0.0123592 *  
## factor(state)25  40.614824  13.285314  3.0571 0.0037925 ** 
## factor(state)26  12.623735  16.515768  0.7643 0.4487419    
## factor(state)27  32.406933  11.771131  2.7531 0.0085468 ** 
## factor(state)28  36.038340  23.184820  1.5544 0.1272547    
## factor(state)29  -8.364926  10.713086 -0.7808 0.4390928    
## factor(state)30 -15.747500  11.752205 -1.3400 0.1871372    
## factor(state)31  27.980646  14.042882  1.9925 0.0525409 .  
## factor(state)32   0.085949  13.169517  0.0065 0.9948223    
## factor(state)33   4.569171  12.207951  0.3743 0.7099950    
## factor(state)34  13.853360   7.413582  1.8686 0.0683386 .  
## factor(state)35  33.850794  10.851760  3.1194 0.0031940 ** 
## factor(state)36  36.923183  15.688738  2.3535 0.0231309 *  
## factor(state)37  71.438932  25.276934  2.8262 0.0070578 ** 
## factor(state)38  -4.283484   8.087885 -0.5296 0.5990396    
## factor(state)39  32.011106  19.872839  1.6108 0.1143756    
## factor(state)40   8.713881   9.190813  0.9481 0.3482530    
## factor(state)41  67.383227  21.355736  3.1553 0.0028907 ** 
## factor(state)42  70.764081  39.542626  1.7896 0.0804080 .  
## factor(state)43  44.381330  11.670909  3.8027 0.0004374 ***
## factor(state)44   1.722018  14.508908  0.1187 0.9060636    
## factor(state)45 -17.753090  14.294312 -1.2420 0.2208275    
## factor(state)46  -3.277062  12.401873 -0.2642 0.7928294    
## factor(year)87   15.402193   7.530318  2.0454 0.0468277 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Odhad pro jednotlivé faktory predstavuje nepozorovaný efekt dané prurezové jednotky (zde ,mesta). Bohuzel tento odhad není konzistentní, viz prednáska. Co je vsak dulezité, odhad parametru pro vliv nezamestnanosti je pri splnení predpokladu pro LSDV estimátor konzistentní. Zároven vidíme, ze je daný parametr i statisticky významný. Tedy s rustem nezamestnanosti o 1 procentní bod vzorste kriminalita o 2,21 zlozinu na 100 tis. obyvatel.

Within estimator

Jelikoz nám LSDV estimátor napáchá trochu neporádek, také prijdeme o znacné mnozství stupnu volnosti, vyuzijeme jiný estimátor FE modelu, který nám bude dávat stejné výsledky. Tímto estimátorem je Fixed Effect estimátor, nebo také within estimátor.

\(crmrte_{it}-\bar{crmrte_i}=\beta_1(unemp_{it}-\bar{unemp}_i)+\gamma d87+c_i-\bar{c}_i+\epsilon_{it}-\bar{\epsilon}_i\)

kdy

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

fe=plm(crmrte~unemp+factor(year),data = crime_panel, model = "within")

bptest(fe)

## 
##  studentized Breusch-Pagan test
## 
## data:  fe
## BP = 6.6535, df = 2, p-value = 0.03591

Jelikoz zamítáme H0 o homoskedasticite pro \(\alpha=0,05\), vyuzijeme robustní odhad pro Var(b). Máme více mozností

# existuje vice moznosti jak odhadnout kovariancni matici Var(b)
coeftest(fe, vcov.=function(x) vcovHC(fe, method="arellano", type="HC1",cluster="group"))    # robustni jak k autokorelaci, tak heteroskedasticite

## 
## t test of coefficients:
## 
##                Estimate Std. Error t value Pr(>|t|)   
## unemp           2.21800    0.80639  2.7505 0.008604 **
## factor(year)87 15.40219    5.12104  3.0076 0.004342 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

coeftest(fe, vcov.=function(x) vcovHC(fe, method="white1", type="HC1"))   # robustni k heteroskedasticite

## 
## t test of coefficients:
## 
##                Estimate Std. Error t value  Pr(>|t|)    
## unemp           2.21800    0.57021  3.8898 0.0003354 ***
## factor(year)87 15.40219    3.62112  4.2534 0.0001082 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Vidíme, ze dostaneme stejný odhad vlivu nezamestnanosti na kriminalitu i casového období.

Nyní nezahrneme casový prvek

fe2=plm(crmrte~unemp,data = crime_panel, model = "within")

bptest(fe2)

## 
##  studentized Breusch-Pagan test
## 
## data:  fe2
## BP = 3.1891, df = 1, p-value = 0.07413

coeftest(fe2, vcov.=function(x) vcovHC(fe, method="white1", type="HC1"))

## 
## t test of coefficients:
## 
##        Estimate Std. Error t value Pr(>|t|)
## unemp -0.018096   0.570206 -0.0317   0.9748

Vidíte, jak je dulezitá specifikace modelu!

First diferences estimator

Bohuzel balícek plm byl predelán a v soucasné dobe chybne vyhodí kazdou promennou, která se chová jako konstanta. To je i prípad casové dummy promenné. Takze nelze zde pouzít FD estimátor.

fd=plm(crmrte~unemp+d87,data = crime_panel, model = "fd")

bptest(fd)

## 
##  studentized Breusch-Pagan test
## 
## data:  fd
## BP = 6.6535, df = 2, p-value = 0.03591

coeftest(fd,vcovHC)

## 
## t test of coefficients:
## 
##        Estimate Std. Error t value Pr(>|t|)
## unemp -0.018096   0.437135 -0.0414   0.9672

fd2=plm(crmrte~unemp,data = crime_panel, model = "fd")

bptest(fd2)

## 
##  studentized Breusch-Pagan test
## 
## data:  fd2
## BP = 3.1891, df = 1, p-value = 0.07413

coeftest(fd2,vcovHC)

## 
## t test of coefficients:
## 
##        Estimate Std. Error t value Pr(>|t|)
## unemp -0.018096   0.437135 -0.0414   0.9672

Random effect model

Opet vyjdeme z modelu

\(crmrte_{it}=\beta_0+\gamma d_{87}+\beta_1 unemp_{it}+c_i+\epsilon_{it}\)

Nyní vsak jiz predpokládáme \(E(unemp_{it}|c_i)=0\). Mohli bychom tak odhadnout rovnici pomocí OLS, jenze odhad by nebyl vydatný. Problémem je následující korelace:

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

pak

\(E(v_{it}, v_{it-1})=E((c_i+\epsilon_{it}),(c_i+\epsilon_{it-1}))=Var(c_i)=\sigma^2_c\)

Odhad naseho RE modelu tak provedeme pomocí metody GLS, respektive FGLS, jelikoz budeme muset odhadnout \(\sigma^2_c\).

\(crmrte_{it}-\lambda \bar{crmrte}_i =\beta_0 + \beta_1 (unemp_{it}-\lambda \bar{unemp}_i)+\gamma d87+(v_{it}-\lambda \bar{v_i})\)

kdy

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

Na upravenou regresní rovnici aplikujeme OLS a toto je pak \(\textbf{RANDOM EFFECT ESTIMATOR}\)

re=plm(crmrte~unemp+factor(year),data = crime_panel, model = "random")

bptest(re)

## 
##  studentized Breusch-Pagan test
## 
## data:  re
## BP = 6.6535, df = 2, p-value = 0.03591

Stále zamítáme H0 o homoskedasticite pro \(\alpha=0,05\). Robustní odhad je napríklad získán jako

coeftest(re, vcovHC(re, method = "arellano", cluester = "group"))

## 
## t test of coefficients:
## 
##                Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)    80.12161    6.51884 12.2908 < 2.2e-16 ***
## unemp           1.74922    0.61813  2.8299  0.005755 ** 
## factor(year)87 13.44965    4.32754  3.1079  0.002529 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Tedy s kazdým dalsím nárustem nezamestnanosti o 1 procentní bod, v prumeru zvýsí kriminalitu o 1,7 cinu na 100 tis. obyvatel. Zároven behem období 82 az 87 vzrostla kriminalita v prumeru o 13,4 cinu na 100 tis. obyvatel.

Vyhodnocení

Máme tedy pracovat spíse s predpoklady FE modelu a nebo RE modelu. Na to nám pomuze najít odpoved Hausmanuv test.

V prípade, ze \(Cor(v_{it},x_{it})= 0\) budou odhady pomocí FE estimátoru, respektive RE estimátoru 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.

Zároven máme predpoklad, ze nebude prítomna heteroskedasticita, ani autokorelace.

Metoda	\(H_0:Exogenní\)	\(H_1:Endogenní\)
FE	Konzistetní a nevydatný	Konzistetní
RE	Konzistetní a vydatný	Nekonzistetní

phtest(fe, re)   # predpoklad homoskedasticity a NEautokorelace

## 
##  Hausman Test
## 
## data:  crmrte ~ unemp + factor(year)
## chisq = 1.9464, df = 2, p-value = 0.3779
## alternative hypothesis: one model is inconsistent

My vsak víme, ze nemuzeme spoléhat na tyto výsledky, jelikoz jsme zamítli hypotézu o homoskedasticite. Autokorelaci netestujeme, jelikoz máme prílis krátký panel T=2

Mozná robustní verze Hausmanova testu

forma=crmrte~unemp+d87                               # bohuzel prikaz faktor zde nefungoval 
phtest(forma, data = crime_panel, method = "chisq",vcov = function(x) vcovHC(x, method="arellano"))

## 
##  Hausman Test
## 
## data:  forma
## chisq = 1.9464, df = 2, p-value = 0.3779
## alternative hypothesis: one model is inconsistent

p-value = 0.3779, tedy nemuzeme zamítnou H0. Z toho muzeme usuzovat, ze zde není problém endogenity a prkloníme se k interpretaci modelu RE.

Priklad 2

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=pdata.frame(kravy, index = c("FARM", "YEAR"))
attach(kravy_panel)

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

Farm = observation, 1-247,
Year = Year, 93-98,
Cows = number of cows,
Land = Amount of land in hectares
Milk = Output in liters
Labor = Number of workers
Feed = Amount of feed used

\(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}=\beta_0+\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

coeftest(lsdv, vcov = parzenHAC)

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

Jelikoz zamítáme H0 pro \(\alpha=0,01\) tedy s velkou pravdepodobností predpokládáme výskyt autokorelace, jiz nemusíme testovat náhodnou slouzku na heteroskedasticitu. Je to z toho duvodu, ze testy na hetersokedastictiu nejsou robustní vuci autokorelaci. Jinak receno v prípade autokorelace nám mohou dávat zavádející výsledky.

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

A provedeme robustní odhad Var(b)

coeftest(fix_wit, vcov.=function(x) vcovHC(x , method="arellano", type="HC1",cluster="group"))

## 
## t test of coefficients:
## 
##          Estimate  Std. Error t value  Pr(>|t|)    
## COWS   3.2686e-02  2.9413e-03 11.1128 < 2.2e-16 ***
## LAND   1.7403e-03  3.4696e-03  0.5016    0.6161    
## LABOR -1.4806e-02  3.5962e-02 -0.4117    0.6806    
## FEED   2.8306e-06  4.0608e-07  6.9706 5.133e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

coeftest(fix_fd, vcov.=function(x) vcovHC(fix_fd , method="arellano", type="HC1",cluster="group"))

## 
## t test of coefficients:
## 
##         Estimate Std. Error t value  Pr(>|t|)    
## COWS  2.8537e-02 2.7162e-03 10.5061 < 2.2e-16 ***
## LAND  3.9670e-03 1.5666e-03  2.5323   0.01146 *  
## LABOR 9.8000e-03 2.6330e-02  0.3722   0.70981    
## FEED  2.3768e-06 3.5015e-07  6.7880 1.763e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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_c}{\sigma^2_c+\sigma^2_{v}})\)

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

kdy

\(\lambda=\dfrac{1}{1+T(\sigma^2_c/\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

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

coeftest(rem, vcov.=function(x) vcovHC(fix_fd , method="arellano", type="HC1",cluster="group"))

## 
## t test of coefficients:
## 
##         Estimate Std. Error t value  Pr(>|t|)    
## COWS  3.5215e-02 2.7162e-03 12.9646 < 2.2e-16 ***
## LAND  1.1965e-03 1.5666e-03  0.7638    0.4451    
## LABOR 3.8866e-02 2.6330e-02  1.4761    0.1401    
## FEED  2.7300e-06 3.5015e-07  7.7965 1.193e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Na rozdíl od FE modelu jiz máme vsechny odhady parametru kladné. To se zdá z ekonomického hlediska logické.

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

Vzdy zacínáme Hausmannovo testem. Jelikoz ten nám pomuze odhalit prípadný problém endogenity. Pokud se nám nepodarí prokázat prítomnost endogenity, pokracujeme dalsími testy.

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

form=log(MILK)~COWS+LAND+LABOR+FEED
phtest(form, data = kravy_panel, method = "aux",vcov = function(x) vcovHC(x, method="arellano", type="HC3"))

## 
##  Regression-based Hausman test, vcov: function(x) vcovHC(x, method
##  = "arellano", type = "HC3")
## 
## data:  form
## chisq = 16.642, df = 4, p-value = 0.002268
## alternative hypothesis: one model is inconsistent

Na hladine významnosti 1%%% zamítáme H0 a kloníme se k pouzití FE modelu a nekterého z jeho robustních estimátoru.

Podle eko. teorie

V minulém príkladu jsme udelali vse jak má být, z ekonometrického pohledu. Problém je v tom, ze ekonomicky nám výsledek nedává smysl. Proto vyjdeme z ekonomické teorie o produkcních funcích a pouzijeme rovnou tu nejznámejsí Cobb-Douglasovu

\(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

pbgtest(Lfix_wit, order = 1)

## 
##  Breusch-Godfrey/Wooldridge test for serial correlation in panel
##  models
## 
## data:  log(MILK) ~ log(COWS) + log(LAND) + log(LABOR) + log(FEED)
## chisq = 0.81881, df = 1, p-value = 0.3655
## alternative hypothesis: serial correlation in idiosyncratic errors

bptest(Lfix_wit)

## 
##  studentized Breusch-Pagan test
## 
## data:  Lfix_wit
## BP = 46.396, df = 4, p-value = 2.037e-09

Autokorelace se nám nepodarila prokázat, avsak H0 o homoskedasticite zamítáme pro \(\alfa=0,01\)

Vyuzijeme nekterý robustní odhad. Zároven se podívejte, jak se mení výsledky.

coeftest(Lfix_wit, vcovSCC(Lfix_wit, type="HC1", maxlag=2))

## 
## t test of coefficients:
## 
##             Estimate Std. Error  t value  Pr(>|t|)    
## log(COWS)  0.6620010  0.0049102 134.8221 < 2.2e-16 ***
## log(LAND)  0.0373525  0.0064900   5.7554  1.09e-08 ***
## log(LABOR) 0.0303995  0.0081175   3.7450 0.0001888 ***
## log(FEED)  0.3825104  0.0237173  16.1279 < 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

Lrem = plm(log(MILK)~log(COWS)+log(LAND)+log(LABOR)+log(FEED), data=kravy_panel, model="random")
pbgtest(Lrem, order = 1)

## 
##  Breusch-Godfrey/Wooldridge test for serial correlation in panel
##  models
## 
## data:  log(MILK) ~ log(COWS) + log(LAND) + log(LABOR) + log(FEED)
## chisq = 40.4, df = 1, p-value = 2.069e-10
## alternative hypothesis: serial correlation in idiosyncratic errors

coeftest(Lrem, vcov.=function(x) vcovHC(Lrem , method="arellano", type="HC1",cluster="group"))

## 
## t test of coefficients:
## 
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.281435   0.110634 47.7379  < 2e-16 ***
## log(COWS)   0.650272   0.027760 23.4252  < 2e-16 ***
## log(LAND)   0.030049   0.013843  2.1707  0.03011 *  
## log(LABOR)  0.035070   0.017324  2.0244  0.04311 *  
## log(FEED)   0.399528   0.015446 25.8666  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A vyhodnotíme podle Hausmanova testu

form=log(MILK) ~ log(COWS)+log(LAND)+log(LABOR)+log(FEED)
phtest(form, data = kravy_panel, method = "chisq" ,vcov = function(x) vcovHC(x, method="arellano", type="HC3"))

## 
##  Hausman Test
## 
## data:  form
## chisq = 12.721, df = 4, p-value = 0.01273
## alternative hypothesis: one model is inconsistent

Pro \(\alpha=0,05\) zamítáme H0 a prikláníme se k modelu FE.

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é.

Pokud nebudete umet casové rady, urcite se nepoustejte do makroekonomických panelu

\(H_0:\) casová rada není stacionární

cipstest(x, lags = 2, type = c("trend", "drift", "none"),
model = c("cmg", "mg", "dmg"),

V prípade existence jednotkového korene, je nutné data tzv. stacionarizovat. Stacionarizaci provedeme pomocí první diference. Opet provedeme test. Pokud jiz budeme moci predpokládat stacionaritu pokracujeme dále

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")) 
coeftest(fe, vcov.=function(x) vcovHC(fe, method="arellano", type="HC1",cluster="group"))

Heteroskedasticity: Random effects

coeftest(rem) 
coeftest(rem, vcovHC) # Heteroskedasticity consistent standard errors
coeftest(rem, vcovHC(rem, method = "arellano", cluester = "group")) 
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é.

Panelova data

Lukas Fryd

3 May 2017

Panelová data

Priklad 1

LSDV

Within estimator

First diferences estimator

Random effect model

Vyhodnocení

Priklad 2

OLS

Fixed effects model

Fixed effects odhadnutý pomocí Least squares dummy variable model

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

Fixed effects - “first difference” FD

Random effects

Který model pouzít?

Fixed model vs. Random model

Podle eko. teorie

Testování modelu na splnení GM

Testování sériové korelace

Testování heteroskedasticity

Testování cross-sectional dependency

Testování stacionarity

Co delat v prípade autokorelace a heteroskedasticity

Heteroskedasticity: Fixed effects

Heteroskedasticity: Random effects

Dodatecný príklad

Priklad

Odhad FE modelu