Summeopgave2192015

## Løsningsforslag til dele af summeopgave fra d. 21.9.2015
require(foreign)

## Loading required package: foreign

## Warning: package 'foreign' was built under R version 3.2.2

require(ggplot2)

## Loading required package: ggplot2

setwd("C:/Users/cmd/Google Drev/RegressionsAnalyse2015/DataVerbeek/ch2")

#1 Indlaeser data fra en fil i STATA format (dvs en fil med extension .dta)
wages <- read.dta(file ="wages1.dta" ) # Husk a R kun søger efter denne fil i den sti der er angivet i setwd
head(wages)

##   exper male school     wage
## 1     9    0     13 6.315296
## 2    12    0     12 5.479770
## 3    11    0     11 3.642170
## 4     9    0     14 4.593337
## 5     8    0     14 2.418158
## 6     9    0     14 2.094058

#2 Laver diverse krydsplot
ggplot(wages,aes(exper,wage))+geom_point()

ggplot(wages,aes(male,wage))+geom_point()

ggplot(wages,aes(school,wage))+geom_point()

# Plot ser ud til at indikere at der et et positivt afkast ved skolegang

#3+4 Det kunne være interessant at estimere afkastet af at tage et års ekstra skolegang. 
# Dette afkast er målt ved beta2. Lad os estimere beta1 og beta2
y    = wages$wage
x    = wages$school

ybar = mean(y)
xbar = mean(x)

ytilde = y-ybar
xtilde = x-xbar

b = sum(xtilde*ytilde)/sum(xtilde^2)
a = ybar-b*xbar

cat(c(a=a,b=b))

## -0.7225066 0.5571617

#5 Tjek af om foc og soc holder
# foc:
foc.a = -2*sum(y-a-b*x)
foc.b = -2*sum(x*(y-a-b*x))
cat(c(foc.a=foc.a,foc.b=foc.b))

## 7.585044e-13 1.055156e-11

# soc
capN   = NROW(y)
soc.aa = 2*capN
soc.ab = 2*sum(x)
soc.ba = 2*sum(x)
soc.bb = 2*sum(x^2)   
soc.hessian = matrix(c(soc.aa,soc.ab,soc.ba,soc.bb),2,2)
# Vi tjekker om soc.hessian er positiv definite ved at beregne eigenvalues som alle skal være positive
er.soc.hessian.pd=eigen(soc.hessian)
cat("Eigenvalues er :",er.soc.hessian.pd$values)

## Eigenvalues er : 915707.8 130.1814

# Da eigenvaerdierne er alle positve kan vi konkludere at soc er opfyldt

#8 Indtegner "regressionslinjen"
## Replicating Figure 2.1 type page 10 in Verbeek
ggplot(wages,aes(school,wage))+geom_point(col="red")+geom_abline(intercept=a,b,col="blue")

#9 Residualplots
# Beregner residualer
e       = y-a-b*x
wages$e = e
wages$i = seq(1,capN)
# Diverse residualplots
ggplot(wages,aes(i,e))+geom_point(col="red")+geom_abline(intercept=0,slope=0,col="blue")

ggplot(wages,aes(school,e^2))+geom_point(col="red")+geom_abline(intercept=0,slope=0,col="blue")

ggplot(wages,aes(exper,e))+geom_point(col="red")+geom_abline(intercept=0,slope=0,col="blue")

ggplot(wages,aes(exper,e^2))+geom_point(col="red")+geom_abline(intercept=0,slope=0,col="blue")

# Det er ikke klart om de estimerede resdualer opfører sig som iid?..eller?

#10 Beregner determinationskoefficienten (R^2)
SST  = sum(ytilde^2)
SSE  = b^2*sum(xtilde^2)
SSR  = sum(e^2)
R.sq = 1-(SSR/SST)
cat("Determinationskoefficienten :", R.sq)

## Determinationskoefficienten : 0.07980201