datam = t(matrix(scan("C:/HY-Data/JPLUOTO/OPETUS/EKONOMETRIAN_PERUSKURSSI/LUENNOT/DATA/cig.txt"),nrow=8))
colnames(datam) = c("year","cpi","pop","packpc","income","tax","avgprs","taxs")


#Ladataan tarvittavat paketit
library(car)
library(sem)


#Nimetään muuttujat

ln_Q_cigar = log(datam[49:96,"packpc"])
ln_P_cigar = log(datam[49:96,"avgprs"]/datam[49:96,"cpi"])
Sales_Tax = datam[49:96,"taxs"]/datam[49:96,"cpi"]
Cig_Tax = datam[49:96,"tax"]/datam[49:96,"cpi"]
ln_Inc = log(datam[49:96,"income"]/(datam[49:96,"pop"]*datam[49:96,"cpi"]))




#1. regression (huom! tulokset eivät vastaa kirjan tuloksia)
model = lm(ln_P_cigar~Sales_Tax)
summary(model)


#S&W 12.10 (huom! tulokset eivät vastaa kirjan tuloksia)
model = tsls(ln_Q_cigar~ln_P_cigar,~Sales_Tax)
summary(model)
#Vert.
Xhat = c(lm(ln_P_cigar~Sales_Tax)$fitted.values)
summary(lm(ln_Q_cigar~Xhat))


#Lisätään kontrolli muuttuja ln_Inc
#S&W 12.15 (huom! tulokset eivät vastaa kirjan tuloksia)
model = tsls(ln_Q_cigar~ln_P_cigar+ln_Inc,~Sales_Tax+ln_Inc)
summary(model)
#Vert.
Xhat = c(lm(ln_P_cigar~Sales_Tax+ln_Inc)$fitted.values)
summary(lm(ln_Q_cigar~Xhat+ln_Inc))


#Lisätään toinen instrumenttimuuttuja Cig_Tax
#S&W 12.16
model = tsls(ln_Q_cigar~ln_P_cigar+ln_Inc,~Cig_Tax+Sales_Tax+ln_Inc)
summary(model)
u_TSLS = model$residuals
#Vert.
#Xhat = c(lm(ln_P_cigar~Cig_Tax+Sales_Tax+ln_Inc)$fitted.values)
#summary(lm(ln_Q_cigar~Xhat+ln_Inc))




#F-testi heikoille instrumenteille
fmodel = lm(ln_P_cigar~Cig_Tax+Sales_Tax+ln_Inc)
#H0: pi_Cig_Tax = pi_Sales_Tax = 0 (vert. S&W yhtälö 12.14)
m = 2                             #(vert S&W yhtälö 12.12)
r = 1                             #(vert S&W yhtälö 12.12)
rm = m + r
q = 2                             #testattavien rajoitteiden lkm.
r = rep(0,q)
R = matrix(0,nc=rm+1,nr=q)
R[1,2] = 1
R[2,3] = 1
linearHypothesis(fmodel,R,r,white.adjust="hc1")




#J-testi 
jmodel = lm(u_TSLS ~ Cig_Tax+Sales_Tax+ln_Inc)
#H0: delta_Cig_Tax = delta_Sales_Tax = 0 (vert. S&W yhtälö 12.17)
m = 2  # (vert S&W yhtälö 12.12)
k = 1  # (vert S&W yhtälö 12.12) 
r = 1
rm = m + r
q = 2
r = rep(0,q)
R = matrix(0,nc=rm+1,nr=q)
R[1,2] = 1
R[2,3] = 1
linearHypothesis(jmodel,R,r,white.adjust="FALSE")

J = m*linearHypothesis(jmodel,R,r,white.adjust="FALSE")$F[2]
pchisq(J,df=m-k,lower.tail = FALSE)






#IV-regressio 10-vuoden-muutoksille

dln_Q_cigar = log(datam[49:96,"packpc"]) - log(datam[1:48,"packpc"])
dln_P_cigar = log(datam[49:96,"avgprs"]/datam[49:96,"cpi"]) - log(datam[1:48,"avgprs"]/datam[1:48,"cpi"])
dSales_Tax = datam[49:96,"taxs"]/datam[49:96,"cpi"] - datam[1:48,"taxs"]/datam[1:48,"cpi"]
dCig_Tax = datam[49:96,"tax"]/datam[49:96,"cpi"] - datam[1:48,"tax"]/datam[1:48,"cpi"]
dln_Inc = log(datam[49:96,"income"]/(datam[49:96,"pop"]*datam[49:96,"cpi"]))-log(datam[1:48,"income"]/(datam[1:48,"pop"]*datam[1:48,"cpi"]))



summary(lm(dln_P_cigar~dln_Inc+dSales_Tax))
F = 12.005^2
model = tsls(dln_Q_cigar~dln_P_cigar+dln_Inc,~dSales_Tax+dln_Inc)
summary(model)


summary(lm(dln_P_cigar~dln_Inc+dCig_Tax))
F = 9.670^2
model = tsls(dln_Q_cigar~dln_P_cigar+dln_Inc,~dCig_Tax+dln_Inc)
summary(model)



model = tsls(dln_Q_cigar~dln_P_cigar+dln_Inc,~dCig_Tax+dSales_Tax+dln_Inc)
u_TSLS = model$residuals
summary(model)

#F-testi heikoille instrumenteille
fmodel = lm(dln_P_cigar~dCig_Tax+dSales_Tax+dln_Inc)
#H0: pi_dCig_Tax = pi_dSales_Tax = 0
m = 2                             
r = 1                             
rm = m + r
q = 2                             
r = rep(0,q)
R = matrix(0,nc=rm+1,nr=q)
R[1,2] = 1
R[2,3] = 1
linearHypothesis(fmodel,R,r,white.adjust="hc1")



#J-testi 
jmodel = lm(u_TSLS ~ dCig_Tax+dSales_Tax+dln_Inc)
#H0: delta_Cig_Tax = delta_Sales_Tax = 0 (vert. S&W yhtälö 12.17)
m = 2
k = 1 
r = 1
rm = m + r
q = 2
r = rep(0,q)
R = matrix(0,nc=rm+1,nr=q)
R[1,2] = 1
R[2,3] = 1
linearHypothesis(jmodel,R,r,white.adjust="FALSE")

J = m*linearHypothesis(jmodel,R,r,white.adjust="FALSE")$F[2]
pchisq(J,df=m-k,lower.tail = FALSE)