#===============================================================================================#
# Metody Ekonometryczne
# Testowanie prostych hipotez
# Jakub Muck 
#===============================================================================================#

#	Wczytanie danych
rm(list=ls())
# Poniżej kod umożliwiający instalacje pakietu AER
# install.packages("AER")
library(AER)
# skorzystaj z wbudowanego zbioru danych
# umożliwiających przeprowadzenie tzw. "growth regressions"
# por. Mankiw, N.G, Romer, D., and Weil, D.N. (1992). A Contribution to the Empirics of Economic Growth. Quarterly Journal of Economics, 107, 407–437.
data(GrowthDJ)

# informacja o danych
help("GrowthDJ")

# Regresja wzrostu bazująca na modelu Solowa
model	=lm(gdpgrowth~log(gdp60)+log(invest)+popgrowth,data=GrowthDJ)
summary(model)

# zastanów się:
# co ta regresja oznacza?
# jaka jest interpretacja parametrów?
# dlaczego gdp60 jest z początku okresu, a zmienna y, czyli gdpgrowth jest z okresu rozpoczynającego się w 1960?
# jaka jest różnica w tempach wzrostu krajów różniących się stopą inwestycji o 10% (czyli np 20% i 22%), ceteris paribus?
# A jak powinna wyglądać specyfikacja gdy chcemy zmierzyć efekt różnicy stopy inwestycji w punktach procentowych?
# co to jest konwergencja i co oznacza spełnienie warunków konwergencji dla znaku parametru przy log(gdp60), czyli
# PKB z początku okresu (wzrost gdpgrowth jest dla okresu 1960-85)

# Liczba stopni swobody
df	=model$df.residual
# Macierz oszacowań (punktowych) oraz błędów standardowych
COEF	=summary(model)$coefficients[,1:2]
COEF

# Rozszerzenie macierzy COEF o:
# (i) oszacowania przedziałowe dla 95% prawdopodobieństwa (1-alpha)
# (ii) statystyki testu t (na istotność oszacowań, tj. H0:beta_i=0)
# (iii) korenspondujace z (ii) p-value
# porównaj Prezentacja 2, slajd 21 
alpha	=0.05
COEF	=cbind(COEF,
            COEF[,1]-COEF[,2]*qt((1-alpha/2),df),
            COEF[,1]+COEF[,2]*qt((1-alpha/2),df),
            COEF[,1]/COEF[,2],
            2*pt(abs(COEF[,1]/COEF[,2]),df,lower.tail=FALSE)
)
colnames(COEF)[3:6]=c("lower bound","upper bound","t","p-value(t)")
COEF

# Prosty test dla jednego współczynnika
# Przykładowo H_0: beta_log(gdp60)=-1
# i H_1: beta_log(gdp60) rożne od -1

statystyka_t=(COEF[2,1]-(-1))/COEF[2,2]
statystyka_t
2*pt(abs(statystyka_t),df,lower.tail=FALSE)
# zastanów się:
# co to jest funkcja pt(.)?
# dlaczego we wzorze jest abs(.)?
# co oznacza lower.tail=FALSE? Co by się działo gdyby parametr ten był ustawiony na FALSE?
# dlaczego w powyższym wzorze jest 2*...?


# Wykorzystanie zmiennych indykatorowych/przekształceń przy pomocy "I()"
# jeśli chcemy stworzyć (tylko lokalnie, na potrzeby regresji, bez zmieniania danych) 
# nową zmienną przy pomocy równania, to możemy wykorzystać
# funkcję I() w formule definiującej model, np dodanie kwadratu to: + I(x^2) +
# I(oecd=="yes") tworzy zmienną binarną dla krajów OECD
model2=lm(gdpgrowth~log(gdp60)+log(invest)+popgrowth+I(oecd=="yes"),data=GrowthDJ)
summary(model2)

# W jaki sposób można oszacować model tylko dla krajów OECD? 
# Są co najmniej 2 sposoby, wskazówka dla jednego z nich - funkcja lm ma opcję subset która powinna wskazywać na wektor 
# binarny wskazujący które obserwacje są włączane do estymacji, np. subset = oecd=="yes" 
# A jak oszacować model kiedy chcemy oszacować model na pełnej próbie, ale zróżnicować parametr związany z konwergencją (i stałą) 
# w grupie OECD-NonOECD?