# 계량경제학 강의 R 코드

# 제8장 다중효과 모형과 그 추정

# 8.4 OLS 추정량이 유일할 조건

# 예제 8.1 여성 더미, 남성 더미, 공선성

data(Wages1,package="Ecdat")
Wages1$male <- as.numeric(Wages1$sex=="male")
Wages1$female <- as.numeric(Wages1$sex=="female")
head(Wages1)
sum(Wages1$male)
sum(Wages1$female)

lm(wage~female,data=Wages1,subset=female==1)  # singular
lm(wage~female,data=Wages1)                   # OK
lm(wage~male+female,data=Wages1)              # singular
lm(wage~male+female-1,data=Wages1)            # OK

# 8.5 더미변수와 상호작용항

data(Wages1, package="Ecdat")
names(Wages1)
Wages1$female <- as.numeric(Wages1$sex=='female')
ols <- lm(log(wage)~school+exper+female,data=Wages1)
summary(ols)

data(Wages1, package="Ecdat")
Wages1$female <- as.numeric(Wages1$sex=='female')
ols <- lm(log(wage)~female*school+exper,data=Wages1)
summary(ols)

data(Wages1, package="Ecdat")
Wages1$female <- as.numeric(Wages1$sex=='female')
eqa <- lm(log(wage)~school,data=Wages1)
eqb <- lm(log(wage)~school+female,data=Wages1)
eqc <- lm(log(wage)~school*female,data=Wages1)
eqd <- lm(log(wage)~school,data=Wages1,subset=female==0)
eqe <- lm(log(wage)~school,data=Wages1,subset=female==1)

# 8.6 제곱항

set.seed(101)
n <- 100
x <- rnorm(n,mean=3.2,sd=1)
u <- rnorm(n,sd=.75)
y <- 3+7*x-x^2+u

plot(y~x,pch=19)

# 예제 8.2 나이의 2차 함수로서의 흡연량

datadir <- "http://econ.korea.ac.kr/~chirokhan/book/data"
smoke <- read.csv(file.path(datadir,"smoke.csv"))
summary(lm(cigs~log(income)+log(cigpric)+educ+age+agesq,data=smoke))

-0.780635/(2*(-0.009106))

summary(smoke$age)

# 예제 8.3 건평과 주택가격

data(Housing,package="Ecdat")
H3 <- Housing[Housing$bedrooms==3,]
ols2 <- lm(log(price)~log(lotsize)+I(log(lotsize)^2),data=H3)
summary(ols2)

-3.29052/(-0.16510*2)

summary(log(H3$lotsize))

ols1 <- lm(log(price)~log(lotsize),data=H3)
summary(ols1)

plot(log(price)~log(lotsize),data=H3)
x <- sort(log(H3$lotsize))
b1 <- ols1$coef
b2 <- ols2$coef
lines(x,b2[1]+b2[2]*x+b2[3]*x^2)
lines(x,b1[1]+b1[2]*x,lty=2)
legend("topleft",c("quadratic","linear"),lty=c(1,2))

# 8.7 맞춘값, 잔차, 제곱합, R제곱

data(Wages1, package="Ecdat")
summary(lm(log(wage)~sex+school+exper,data=Wages1))

# 조정된 R제곱

data(Wages1, package="Ecdat")
Wages1$rnd <- rnorm(nrow(Wages1))
summary(lm(log(wage)~sex+school+exper,data=Wages1))$r.sq
summary(lm(log(wage)~sex+school+exper+rnd,data=Wages1))$r.sq

# 8.A 최소제곱 추정량의 행렬 표현

# 예제 8.4 수동식 OLS 추정

data(Housing,package="Ecdat")
lm(log(price)~log(lotsize)+bedrooms,data=Housing)

Y <- log(Housing$price)
X <- cbind(1,log(Housing$lotsize),Housing$bedrooms)
solve(t(X)%*%X,t(X)%*%Y)