Tietokoneharjoitus 4

              storage  display     value
variable name   type   format      label      variable label
-------------------------------------------------------------------------------
county          int    %9.0g                  county identifier
year            byte   %9.0g                  81 to 87
crmrte          float  %9.0g                  crimes committed per person
prbarr          float  %9.0g                  'probability' of arrest
prbconv         float  %9.0g                  'probability' of conviction
prbpris         float  %9.0g                  'probability' of prison sentenc
avgsen          float  %9.0g                  avg. sentence, days
polpc           float  %9.0g                  police per capita
density         float  %9.0g                  people per sq. mile
taxpc           float  %9.0g                  tax revenue per capita
west            byte   %9.0g                  =1 if in western N.C.
central         byte   %9.0g                  =1 if in central N.C.
urban           byte   %9.0g                  =1 if in SMSA
pctmin80        float  %9.0g                  perc. minority, 1980
wcon            float  %9.0g                  weekly wage, construction
wtuc            float  %9.0g                  wkly wge, trns, util, commun
wtrd            float  %9.0g                  wkly wge, whlesle, retail trade
wfir            float  %9.0g                  wkly wge, fin, ins, real est
wser            float  %9.0g                  wkly wge, service industry
wmfg            float  %9.0g                  wkly wge, manufacturing
wfed            float  %9.0g                  wkly wge, fed employees
wsta            float  %9.0g                  wkly wge, state employees
wloc            float  %9.0g                  wkly wge, local gov emps
mix             float  %9.0g                  offense mix: face-to-face/other
pctymle         float  %9.0g                  percent young male
d82             byte   %9.0g                  =1 if year == 82
d83             byte   %9.0g                  =1 if year == 83
d84             byte   %9.0g                  =1 if year == 84
d85             byte   %9.0g                  =1 if year == 85
d86             byte   %9.0g                  =1 if year == 86
d87             byte   %9.0g                  =1 if year == 87
lcrmrte         float  %9.0g                  log(crmrte)
lprbarr         float  %9.0g                  log(prbarr)
lprbconv        float  %9.0g                  log(prbconv)
lprbpris        float  %9.0g                  log(prbpris)
lavgsen         float  %9.0g                  log(avgsen)
lpolpc          float  %9.0g                  log(polpc)
ldensity        float  %9.0g                  log(density)
ltaxpc          float  %9.0g                  log(taxpc)
lwcon           float  %9.0g                  log(wcon)
lwtuc           float  %9.0g                  log(wtuc)
lwtrd           float  %9.0g                  log(wtrd)
lwfir           float  %9.0g                  log(wfir)
lwser           float  %9.0g                  log(wser)
lwmfg           float  %9.0g                  log(wmfg)
lwfed           float  %9.0g                  log(wfed)
lwsta           float  %9.0g                  log(wsta)
lwloc           float  %9.0g                  log(wloc)
lmix            float  %9.0g                  log(mix)
lpctymle        float  %9.0g                  log(pctymle)
lpctmin         float  %9.0g                  log(pctmin)
clcrmrte        float  %9.0g                  lcrmrte - lcrmrte[_n-1]
clprbarr        float  %9.0g                  lprbarr - lprbarr[_n-1]
clprbcon        float  %9.0g                  lprbconv - lprbconv[_n-1]
clprbpri        float  %9.0g                  lprbpri - lprbpri[t-1]
clavgsen        float  %9.0g                  lavgsen - lavgsen[t-1]
clpolpc         float  %9.0g                  lpolpc - lpolpc[t-1]
cltaxpc         float  %9.0g                  ltaxpc - ltaxpc[t-1]
clmix           float  %9.0g                  lmix - lmix[t-1]

Tehtävä 5

Estimoi lineaarinen malli POLS estimaattorilla (pooled ordinary least squares) kayttaen kaikkia vuosia 81-87. Mallissa vastemuuttuja on log(crmrte) ja selittavat muuttujat ovat log(prbarr), log(prbconv), log(prbpris), log(avgsen) ja log(polpc).

file<-"http://cc.oulu.fi/~jklemela/panel/cornwell.raw"
data<-read.table(file=file)

y<-log(data[,3])

x1<-log(data[,4])
x2<-log(data[,5])
x3<-log(data[,6])
x4<-log(data[,7])
x5<-log(data[,8])

reg.model<-lm(y ~ x1+x2+x3+x4+x5)

summary(reg.model)

Call:
lm(formula = y ~ x1 + x2 + x3 + x4 + x5)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.93632 -0.19952  0.02951  0.25318  1.29941 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -2.20673    0.23869  -9.245  < 2e-16 ***
x1          -0.72151    0.03671 -19.655  < 2e-16 ***
x2          -0.54928    0.02627 -20.909  < 2e-16 ***
x3           0.23797    0.06643   3.582 0.000367 ***
x4          -0.06520    0.05535  -1.178 0.239281    
x5           0.36252    0.02996  12.100  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.3789 on 624 degrees of freedom
Multiple R-squared: 0.5658,	Adjusted R-squared: 0.5624 
F-statistic: 162.7 on 5 and 624 DF,  p-value: < 2.2e-16 

######################################

K<-6
n<-length(y)
x<-matrix(0,n,K)
x[,1]<-1
x[,2]<-x1
x[,3]<-x2
x[,4]<-x3
x[,5]<-x4
x[,6]<-x5

A<-t(x)%*%x
invA<-solve(A,diag(1,K))
b<-invA%*%t(x)%*%y
b

[1,] -2.20673453
[2,] -0.72151130
[3,] -0.54927675
[4,]  0.23797055
[5,] -0.06519926
[6,]  0.36252309

###########################################

county<-data[,1]
year<-data[,2]
uc<-unique(county)
uy<-unique(year)
N<-length(uc)
T<-length(uy)

residu<-matrix(0,N,T)
for (i in 1:N){
    xi<-x[((i-1)*7+1):(i*7),]
    yi<-y[((i-1)*7+1):(i*7)]
    residu[i,]<-yi-xi%*%b
}

Ahat<-matrix(0,K,K)
for (i in 1:N){
    xi<-x[((i-1)*7+1):(i*7),]
    Ahat<-Ahat+t(xi)%*%xi
}
Ahat<-Ahat/N

Bhat<-matrix(0,K,K)
for (i in 1:N){
    xi<-x[((i-1)*7+1):(i*7),]
    ui<-matrix(residu[i,],T,1)
    Bhat<-Bhat+t(xi)%*%ui%*%t(ui)%*%xi
}
Bhat<-Bhat/N

sighat2<-mean(residu^2)
invAhat<-solve(Ahat,diag(1,K))
avar.robust<-invAhat%*%Bhat%*%invAhat/N
avar<-sighat2*invAhat/N

###################################


sdhats<-matrix(0,K,1)
tstats<-matrix(0,K,1)
pvals<-matrix(0,K,1)
for (k in 1:K){
   sdhats[k]<-sqrt(avar[k,k])
   tstats[k]<-b[k]/sdhats[k]
   pvals[k]<-2*(1-pnorm(abs(tstats[k])))
}
sdhats
tstats
pvals

> sdhats
           [,1]
[1,] 0.23755308
[2,] 0.03653369
[3,] 0.02614469
[4,] 0.06611309
[5,] 0.05508736
[6,] 0.02981777
> tstats
           [,1]
[1,]  -9.289438
[2,] -19.749203
[3,] -21.009110
[4,]   3.599447
[5,]  -1.183561
[6,]  12.157956
> pvals
            [,1]
[1,] 0.000000000
[2,] 0.000000000
[3,] 0.000000000
[4,] 0.000318895
[5,] 0.236586851
[6,] 0.000000000

# robustit ####################
sdhats<-matrix(0,K,1)
tstats<-matrix(0,K,1)
pvals<-matrix(0,K,1)
for (k in 1:K){
   sdhats[k]<-sqrt(avar.robust[k,k])
   tstats[k]<-b[k]/sdhats[k]
   pvals[k]<-2*(1-pnorm(abs(tstats[k])))
}
sdhats
tstats
pvals

> sdhats
           [,1]
[1,] 0.85072888
[2,] 0.10845868
[3,] 0.06977112
[4,] 0.10549265
[5,] 0.10197897
[6,] 0.11856455
> tstats
           [,1]
[1,] -2.5939340
[2,] -6.6524070
[3,] -7.8725522
[4,]  2.2558022
[5,] -0.6393402
[6,]  3.0576012
> pvals
             [,1]
[1,] 9.488471e-03
[2,] 2.883382e-11
[3,] 3.552714e-15
[4,] 2.408302e-02
[5,] 5.226016e-01
[6,] 2.231163e-03

#################################

OLS-kaavat

Ahatb<-t(x)%*%x/(N*T)
sighat2<-sum(residu^2)/(N*T-K)
invAhatb<-solve(Ahatb,diag(1,K))
avarb<-sighat2*invAhatb/(N*T)
sdhats<-matrix(0,K,1)
tstats<-matrix(0,K,1)
pvals<-matrix(0,K,1)
for (k in 1:K){
   sdhats[k]<-sqrt(avarb[k,k])
   tstats[k]<-b[k]/sdhats[k]
   pvals[k]<-2*(1-pnorm(abs(tstats[k])))
}
sdhats
tstats
pvals

> sdhats
           [,1]
[1,] 0.23869243
[2,] 0.03670891
[3,] 0.02627009
[4,] 0.06643018
[5,] 0.05535157
[6,] 0.02996078
> tstats
           [,1]
[1,]  -9.245097
[2,] -19.654934
[3,] -20.908828
[4,]   3.582265
[5,]  -1.177912
[6,]  12.099923
> pvals
             [,1]
[1,] 0.0000000000
[2,] 0.0000000000
[3,] 0.0000000000
[4,] 0.0003406274
[5,] 0.2388318499
[6,] 0.0000000000