# Econ 512   HW 2
# Jiehong Qiu
# 10 Nov. 2016


# Q1,a
rm(list = ls())
setwd("C:/Users/qjhws/Desktop/Classes/512TA/HW2")

set.seed(1234)

Ndraw <- 10000
x  <- matrix(0,nrow=Ndraw,ncol=1)
y  <- matrix(0,nrow=Ndraw,ncol=1)
mean_x <- 0.5
mean_y <- -0.5
sigma  <- 1
rho <- -0.15

x_0 <- 0
for (i in 1:Ndraw){
    if(i == 1) {
        y[i] <- rnorm(1,mean_y+rho*solve(sigma)*(x_0-mean_x),
        sqrt(sigma-rho*solve(sigma)*rho))
    }
  
    else{
    j <- i - 1
    y[i] <- rnorm(1,mean_y+rho*solve(sigma)*(x[j]-mean_x),
              sqrt(sigma-rho*solve(sigma)*rho))
    }

    x[i] <- rnorm(1,mean_x+rho*solve(sigma)*(y[i]-mean_y),
                    sqrt(sigma-rho*solve(sigma)*rho))
}
hist(x,freq=FALSE)
hist(y,freq=FALSE)


# Q1,b
set.seed(1234)
df <- 5
chi_RN <- rchisq(Ndraw,df)
hist(chi_RN,freq=FALSE)


# Q1,c

y_c <- rnorm(Ndraw,2,sqrt(0.9))
x_c <- (1/sqrt(0.9))*(y_c-2)
mean(x_c)

df <- 100
z <- x_c^2
i <- 1
while(i<= df){
  z <- z + ((1/sqrt(0.9))*(rnorm(Ndraw,2,sqrt(0.9))-2))^2
  i <- i + 1
}

T <- x_c/(sqrt(z/df))
hist(T,freq=FALSE)


# Q1,d 
set.seed(123)
p <- rpois(Ndraw, 4)
hist(p,freq=FALSE)

# Q1,e
set.seed(123)
MN <- rmultinom(10, size = 1, prob=c(1/5,1/10,1/2,1/5))




#Q2---------------------------------------------------------------------------------
# a
set.seed(123)
x_bar <- rnorm(10,0,sqrt(6))
theta_n <- exp(mean(x_bar))
bias <- theta_n-1
bias


# b
bias_est <- c()  # create a empty vector
Nsample <- 100   # number of sampling from x_bar

i <- 1
while(i <= Nsample){
    mean_bs <- mean(sample(x_bar, 10,replace = TRUE))
    tem_bias <- mean_bs - mean(x_bar)
    bias_est <- rbind(bias_est,tem_bias)
    i <- i + 1
}
bias_est
bias_est <- mean(bias_est)
theta_nStar <- theta_n-bias_est
theta_nStar
bs_bias <- theta_nStar- 1  # bias of the bias-correted estimator
bs_bias


#Q3--------------------------------------------------------------------------------
set.seed(123)
x <- c()
y <- c()
for (i in 1:Ndraw){
  x_temp  <- runif(1, 0, 1)
  y_temp  <- runif(1, 0, 1)
      if(((x_temp-0.5)^2+(y_temp-0.5)^2) <= 0.25 ){
      x <- rbind(x,x_temp)
      y <- rbind(y,y_temp)
      }
  }

result <- sin(sqrt(abs(log(x+y))))
result <- sum(result)/nrow(x)


#Q4--------------------------------------------------------------------------------

data1 <- read.table( "time_var.txt", header = TRUE)
data2 <- read.table( "time_invar.txt", header = TRUE)
data <- merge(data1,data2, by=c("PersonID"))
attach(data)




#Model 1 
x <- cbind(rep(1,nrow(data)),Education,	Experience,	TimeTrend)
y <- LogWage

bols <- solve((t(x)%*%x),(t(x)%*%y))  #ols estimators

n <- nrow(x)
k <- ncol(x)
e <- y-x%*%bols # residual
s2 <- (t(e)%*%e)/(n-k) 
vb <- s2[1,1]*solve(t(x)%*%x)
se <- sqrt(diag(vb))
tval=bols/se

tt <- data.frame(col1=bols,
                 col2=se,
                 col3=tval) 

colnames(tt) <- c("estimate","s.e.","t")
rownames(tt) <- c("constant", "educ", "exper","TimeTrend")
print(tt)


#Model 2

x <- cbind(rep(1,nrow(data)),Education,	Experience,	TimeTrend, 
           Ability,	MotherEducation, FatherEducation,	BrokenHome,	NumberSiblings)

y <- LogWage

bols <- solve((t(x)%*%x),(t(x)%*%y))  #ols estimators

n <- nrow(x)
k <- ncol(x)
e <- y-x%*%bols # residual
s2 <- (t(e)%*%e)/(n-k) 
vb <- s2[1,1]*solve(t(x)%*%x)
se <- sqrt(diag(vb))
tval=bols/se

tt <- data.frame(col1=bols,
                 col2=se,
                 col3=tval) 

colnames(tt) <- c("estimate","s.e.","t")
rownames(tt) <- c("constant", "educ", "exper","TimeTrend","Ability",	
                  "MotherEducation", "FatherEducation",	"BrokenHome",	"NumberSiblings")
print(tt)



#Model 3
x <- cbind(rep(1,nrow(data)),Education,	Experience,	TimeTrend, 
           Ability,	MotherEducation, FatherEducation,	BrokenHome,	NumberSiblings)

y <- LogWage

data <- cbind(PersonID,LogWage,rep(1,nrow(data)),TimeTrend,Experience)
data_FE <- aggregate(data[,2:5],list(PersonID),mean)
colnames(data_FE) <- c("PersonID","LogWage_M","Cons_M","Time_M","Exper_M")
Final_data <- merge(data,data_FE, by=c("PersonID"))
attach(Final_data)
x <- cbind((TimeTrend-Time_M),(Experience-Exper_M))
y <- LogWage-LogWage_M

bols <- solve((t(x)%*%x),(t(x)%*%y))  #ols estimators

n <- nrow(x)
k <- ncol(x)
e <- y-x%*%bols # residual
s2 <- (t(e)%*%e)/(n-k-nrow(data_FE))  # the degree of freedom is nT-n-k. 
vb <- s2[1,1]*solve(t(x)%*%x)
se <- sqrt(diag(vb))
tval=bols/se

tt <- data.frame(col1=bols,
                 col2=se,
                 col3=tval) 

colnames(tt) <- c("estimate","s.e.","t")
rownames(tt) <- c( "educ", "exper")
print(tt)