AKM cvicení 3 Zobecnený lineární model

Gauss-Markov teorém

Lineární v parametrech
\(E(\epsilon)=0\) toto znáte
\(E(\epsilon|X)=E(\epsilon)=0\)
Neexistuje dokonalá multikolinearita
\(Var(\epsilon|X)=\sigma^2 I\)

Po splnení predpokladu získáme odhad, který je \(\textbf{BLUE}\).

Odhad je tak:

Nezkreslený
Konzistentní
Vydatný

Model

Nasimulujeme si populacní model:

\(wage=-50+10IQ+50educ+3expert+\epsilon\) kdy \(\epsilon \sim N(0,1)\)

V následujícím cvicení se bude jednat o základní model, u kterého vsak budeme postupne menit nekteré parametry.

Model 1

Jedná se o model splòující GM. Bude se jednat o benchmark v nasí analýze.

set.seed(100)
sigma=1
N=10000
IQ=rnorm(N,100,5)               # Vygenerujeme hodnoty IQ z rovnomerného rozdslení v rozmezí 0 -30 
expert=runif(N, min=0, max=30)  # Vygenerujeme hodnoty expet z rovnomerného rozdelení v rozmezí 0 -30  
educ = runif(N, min=9, max=20)  # Vygenerujeme hodnoty educ z rovnomerného rozdelení v rozmezí 9-20 
nevim1=rnorm(N,10,50) 
nevim2=rnorm(N,10,50) 
nevim3=rnorm(N,10,50) 
nevim4=rnorm(N,10,50) 
nevim5=rnorm(N,10,50) 
nevim6=rnorm(N,10,50) 
nevim7=rnorm(N,10,50) 
nevim8=rnorm(N,10,50) 
nevim9=rnorm(N,10,50) 
nevim10=rnorm(N,10,50) 
e = rnorm(N, 0, sigma)          # vygenerujeme náhodnou slozku se strední hodnotou 0 a rozptylem 1
wage = -50 +10*IQ+ 50*educ+3*expert + e # Populacní regresní funkce, proces, který generuje velikost wage

populace1 = data.frame(IQ,expert, educ, wage,nevim1,nevim2,nevim3,nevim4,nevim5,nevim6,nevim7,nevim8,nevim9,nevim10)  # Ulozíme data do data.frame

Odhadnete regresní model:

\(wage=\beta_0+\beta_1IQ+\beta_2educ+\beta_3expert+\beta_4nevim1+...+\beta_{13}nevim10+\epsilon\)

regrese=lm(wage~IQ+educ+expert+nevim1+nevim2+nevim3+nevim4+nevim5+nevim6+nevim7+nevim8+nevim9+nevim10, data = populace1)
summary(regrese)

## 
## Call:
## lm(formula = wage ~ IQ + educ + expert + nevim1 + nevim2 + nevim3 + 
##     nevim4 + nevim5 + nevim6 + nevim7 + nevim8 + nevim9 + nevim10, 
##     data = populace1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.6237 -0.6683 -0.0137  0.6731  3.6289 
## 
## Coefficients:
##               Estimate Std. Error   t value Pr(>|t|)    
## (Intercept) -4.988e+01  2.057e-01  -242.516   <2e-16 ***
## IQ           9.999e+00  1.999e-03  5001.499   <2e-16 ***
## educ         5.000e+01  3.135e-03 15950.165   <2e-16 ***
## expert       2.999e+00  1.148e-03  2613.402   <2e-16 ***
## nevim1      -1.547e-04  1.973e-04    -0.784   0.4330    
## nevim2      -1.889e-04  2.003e-04    -0.943   0.3458    
## nevim3       1.712e-04  1.989e-04     0.861   0.3895    
## nevim4       1.220e-04  1.991e-04     0.613   0.5399    
## nevim5       3.822e-04  1.990e-04     1.921   0.0548 .  
## nevim6       6.063e-06  1.987e-04     0.031   0.9757    
## nevim7      -8.877e-05  1.980e-04    -0.448   0.6540    
## nevim8       2.191e-04  2.002e-04     1.095   0.2738    
## nevim9       1.348e-04  2.000e-04     0.674   0.5001    
## nevim10      1.816e-04  2.007e-04     0.905   0.3657    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9939 on 9986 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:      1 
## F-statistic: 2.22e+07 on 13 and 9986 DF,  p-value: < 2.2e-16

co je p-value?
F-test

Omezený vs. neomezený model.

\(H_0: \beta_4=\beta_5=...=\beta_13=0\) vs. \(H_1: non H_0\)

anova(omezeny,neomezeny)

\(F=\dfrac{SSR_O-SSR_{NE}}{SSR_{NE}} \dfrac{n-k-1}{q} \sim F[q,n-k-1]\)

Reálná data

install.packages("AER")  # nainstalovat balícek

library(plm)     # nainstalovaný balícek musíme vyvolat

## Warning: package 'plm' was built under R version 3.2.5

## Loading required package: Formula

data("Gasoline")

head(Gasoline)

##   country year lgaspcar  lincomep      lrpmg  lcarpcap
## 1 AUSTRIA 1960 4.173244 -6.474277 -0.3345476 -9.766840
## 2 AUSTRIA 1961 4.100989 -6.426006 -0.3513276 -9.608622
## 3 AUSTRIA 1962 4.073177 -6.407308 -0.3795177 -9.457257
## 4 AUSTRIA 1963 4.059509 -6.370679 -0.4142514 -9.343155
## 5 AUSTRIA 1964 4.037689 -6.322247 -0.4453354 -9.237739
## 6 AUSTRIA 1965 4.033983 -6.294668 -0.4970607 -9.123903

lgaspcar logarithm of motor gasoline consumption per car
lincomep logarithm of real per-capita income
lrpmg logarithm of real motor gasoline price
lcarpcap logarithm of the stock of cars per capita

regrese2=lm(lgaspcar~lincomep+lrpmg+lcarpcap,data=Gasoline)
summary(regrese2)

## 
## Call:
## lm(formula = lgaspcar ~ lincomep + lrpmg + lcarpcap, data = Gasoline)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.38412 -0.15307 -0.04981  0.16529  0.59684 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.39133    0.11693   20.45   <2e-16 ***
## lincomep     0.88996    0.03581   24.86   <2e-16 ***
## lrpmg       -0.89180    0.03031  -29.42   <2e-16 ***
## lcarpcap    -0.76337    0.01861  -41.02   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.21 on 338 degrees of freedom
## Multiple R-squared:  0.8549, Adjusted R-squared:  0.8536 
## F-statistic:   664 on 3 and 338 DF,  p-value: < 2.2e-16

umet interpretovat odhady
vedet rozdíl mezi \(R^2\) a \(Adj \ R^2\)
co ríká F-test?

Nekdy nás zajímá testování “nejaké” teorie. Zde nás zajímá, zdali je “poptávka” jednotkove elastická. PS poptávka se odhaduje jinak, pozor na to. \(H_0: \beta_2=-1\) vs. \(H_1: \beta_2 \neq 0\)

\(t_j=\dfrac{b_j-\beta_j}{sd(b_j)} \sim t_{n-k-1}\)

length(Gasoline$lgaspcar)
qt()

Porusení nekterých Gauss-Markov predpokladu

V následující cásti si zobrazíme, co se deje s odhady parametru v prípade:

endogenity
heteroskedasticity

Porusení exogenity promenných

Predpoklad exogenity nám ríká \(E(\epsilon|X)=E(\epsilon)=0\)

V prípade jeho porusení, pak platí, ze:

\(E(\epsilon|X) \neq 0\)

Model 2

Nastane situace, kdy bude náhodná slozka \(\epsilon\) korelována s IQ. \(\textbf{POZNAT Z KODU PORUSENÍ EXOGENITY}\)

N=10000
v=rnorm(N,0,1)
w=rnorm(N,0,1)
IQ=100+10*v       # IQ obsahuje slozku v, která je i v náhodné sloZce
educ=runif(N, min=9, max=20) 
expert=runif(N, min=0, max=30)
e2 =0+10*v+10*sqrt(1-0.8^2)*w      # vygenereujeme náhodnou sloZku, která je korelovaná s náhodnou slozkou v promenné IQ 
wage2 = -50 +10*IQ+ 50*educ+3*expert + e2 #
populace2 = data.frame(IQ,expert, educ, wage2)

Odhad parametru je zkreslený a nekonzistetní

Zkreslený/nezkreslený odhad simulace, KOD NEMUSÍTE ZNÁT

Nejrpve odhadneme model 1 R-krát

R = 10^4  # R = pocet replikací; uloZíme si do promenné, aby to slo pozdeji snadno zmenit.

beta.1.NEZ = rep(NA, R)
beta.2.NEZ = rep(NA, R)
beta.3.NEZ = rep(NA, R)

for (r in 1:R) {  # Pro kazdou replikaci r = 1 az R:
  model.1 = lm(wage ~ IQ+educ+expert, data=populace1, subset=sample(N, 50, replace=TRUE))  # Náhodne vybereme 50 hodnot z populace a spocteme regresní model.
# for cyklus zajistí, ze se tento proces bude opakovat R-krát  
    # Ulozíme odhadnuté parametry.
  beta.1.NEZ[r] = model.1$coef[2]
  beta.2.NEZ[r]=model.1$coef[3]
  beta.3.NEZ[r]=model.1$coef[4]
}

Dále odhadneme model 2 R-krát.

R = 10^4  # R = pocet replikací; ulozíme si do promenné, aby to slo pozdeji snadno zmenit.
beta.1.ZK = rep(NA, R)
beta.2.ZK = rep(NA, R)
beta.3.ZK = rep(NA, R)

for (r in 1:R) {  # Pro kazdou replikaci r = 1 az R:
  model.2 = lm(wage2 ~ IQ+educ+expert, data=populace2, subset=sample(N, 50, replace=TRUE))  # Náhodne vybereme 50 hodnot z populace a spocteme regresní model.
# for cyklus zajistí, ze se tento proces bude opakovat R-krát  
    # Ulozíme odhadnuté parametry.
  beta.1.ZK[r] = model.2$coef[2]
  beta.2.ZK[r]=model.2$coef[3]
  beta.3.ZK[r]=model.2$coef[4]
}

Vybereme si nekterý z odhadu, napríklad odhad \(\beta_1\) a odhadneme jeho rozdelení jak pro model \(1\), tak pro model \(2\). Model M1 splòuje GM. Model M2 porusuje predpoklad exogenity.

beta.vydat = c(beta.1.NEZ, beta.1.ZK)  # Vytvoríme vektor odhadu bety
koef.vydat = c(rep("beta.1.NEZ", R),    # pro vytvorení názvu, vektor, ve kterém replikujeme název R-krát
          rep("beta.1.ZK", R))

data.zkresleny = data.frame(beta.vydat, koef.vydat)  # Vytvoríme data frame, spárujeme název s daty

Hodnoty z data.frame si zobrazíme pomocí neparametrického odhadu hustoty pravdepodobnosti.

install.packages("ggplot2")

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 3.2.5

graf = ggplot(data.zkresleny, aes(beta.vydat, fill=koef.vydat, color=koef.vydat)) +  
  geom_density(alpha=.4) +  # Nahradíme histogram neparametrickým odhadem hustoty.
  xlab("odhadovaná hodnota") +
  ylab("hustota")+
  xlim(c(9,12))  # Omezíme rozsah osy x
graf

Konzistetní/nekonzistetní odhad

Simulace konzistentního odhadu, bude nás zajímat pouze odhad \(\beta_1\).

R=5000                # Maximální velikost výberu z populace
beta.1.M3 = rep(NA, R)

for (r in 20:R) {
    model.3 = lm(wage ~ IQ+educ+expert, data=populace1, subset=sample(N, r, replace=TRUE) )
   beta.1.M3[r] = model.3$coef[2]
}

Simulace nekonzistentního odhadu, bude nás zajímat pouze odhad \(\beta_1\).

R=5000                # Maximální velikost výberu z populace
beta.1.M4 = rep(NA, R)

for (r in 20:R) {
    model.4 = lm(wage2 ~ IQ+educ+expert, data=populace2, subset=sample(N, r, replace=TRUE) )
   beta.1.M4[r] = model.4$coef[2]
}

A data si zobrazíme v grafu

par(mfrow=c(2,1))
plot.ts(beta.1.M3,col="red")
plot.ts(beta.1.M4,col="green")

Caption for the picture.

Heteroskedasticita

\(Var(\epsilon|X)=\sigma^2 I\) homoskedasticita

\(Var(\epsilon|X)=\sigma^2 \Omega\) heteroskedasticita

\(Var(\epsilon|X)=\sigma^2 h(X)\)

Pro jednodusí predstavu uvazujme prímkovou regresi

\(Var(\epsilon|x)=\sigma^2 h(x)\) kde \(h(x)\) predstavuje funkci “typ” heteroskedasticity. Rozptyl musí být vzdy nezáporný!

Model heteroskedasticita

\(\textbf{POZNAT Z KODU HETEROSKEDASTICITU}\)

N=10000
IQ=rnorm(N,100,5)
expert=runif(N, min=0, max=30)
educ = runif(N, min=9, max=20)  
sigma=2*log(educ)              # tvar podmíneného rozptylu, rozptyl závisí na jednom z regresoru 
e4= rnorm(N, 0, sigma)  
wage4 = -50 +10*IQ+ 50*educ+3*expert + e4       
populace5 = data.frame(IQ,expert, educ, wage4)

Jen pro zajímavost si odhadneme model.

regrese_hetero=lm(wage4~IQ+educ+expert,data = populace5)   # provedeme odhad regresní rovnice, nyní nedeláme výber
summary(regrese_hetero)

## 
## Call:
## lm(formula = wage4 ~ IQ + educ + expert, data = populace5)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -18.6061  -3.5705  -0.0338   3.5245  21.7284 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -48.949478   1.101469  -44.44   <2e-16 ***
## IQ            9.988045   0.010680  935.21   <2e-16 ***
## educ         50.011133   0.016865 2965.39   <2e-16 ***
## expert        2.997855   0.006179  485.20   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.339 on 9996 degrees of freedom
## Multiple R-squared:  0.999,  Adjusted R-squared:  0.999 
## F-statistic: 3.297e+06 on 3 and 9996 DF,  p-value: < 2.2e-16

Jsou parametry signifikantní?

Test na heteroskedasticitu

Pouzijeme napríklad Breusch-Pagan test. Ten je zalozen na následujím predpokladu.

\(\epsilon^2=\gamma_0+\gamma_1 x_1+\gamma_2 x_2+...+\gamma_k x_k + v\)

Pokud by byl nekterý parametr signifikantní, tak poté muzeme mluvit o heteroskedasticite. Je treba otestovat parametry. Jednou z mozností je napríklad F-test, nebo LM-test.

\(H_0: \gamma_1=\gamma_2=...=\gamma_k=0\)

\(H_1: neplatí \ H_0\)

install.packages("lmtest")

library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

bptest(wage4~IQ+educ+expert)

## 
##  studentized Breusch-Pagan test
## 
## data:  wage4 ~ IQ + educ + expert
## BP = 165.3, df = 3, p-value < 2.2e-16

Nyní si zkusíme sami sestrojit BP test.

reg_hetero=lm(wage4~IQ+educ+expert)
resid_hetero=reg_hetero$residuals          # ulozíme si residua do promenné resid_hetero

reg_res=lm(resid_hetero^2~IQ+educ+expert) # provedeme regresi ve stylu BP testu
summary(reg_res)

## 
## Call:
## lm(formula = resid_hetero^2 ~ IQ + educ + expert)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -38.35 -23.99 -14.85   8.13 434.47 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.00286    8.45673   0.473    0.636    
## IQ           0.01050    0.08200   0.128    0.898    
## educ         1.67325    0.12948  12.922   <2e-16 ***
## expert      -0.04846    0.04744  -1.022    0.307    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 40.99 on 9996 degrees of freedom
## Multiple R-squared:  0.01653,    Adjusted R-squared:  0.01624 
## F-statistic:    56 on 3 and 9996 DF,  p-value: < 2.2e-16

Pokud není splnena podmínka homoskedasticity, metoda OLS jiz neposkytuje BLUE odhad. Jinými slovy, odhad pomocí OLS jiz není vydatný. Existuje jiná metoda, která bude BLUE. Odhad je stále NEZKRESLENÝ A KONZISTENTNÍ, pokud jsou splneny ostatní podmínky!!!

Co delat v prípade prítomnosti heteroskedasticity?