rm(list=ls())
library(foreign)
setwd("/Users/innoshuadmin/Dropbox/RDHeteroTest/RDMonoTest")
romanian <- read.dta(file="./APP1/romanian_twoschooltown.dta")
data<- romanian[c("bcg","bcg_Dmean","bct","bct_Dmean","agus","dzag","zga","meanzag_jk","bct","town","Y","attend","minscore_bytownbyyear_Dmean","maxscore_bytownbyyear_Dmean")]
newdata <-subset(data,is.na(bcg)==0)
set.seed(123)


y=data$bcg_Dmean
x=data$dzag
w=data$meanzag_jk
t=data$attend
w=(w-min(w))/(max(w)-min(w))
wx=w*x
xsq=x^2
wsq=w^2 
pos=(x>0)*1
xpos=x*pos
wpos=w*pos
wxpos=wx*pos
xsqpos=xsq*pos
wsqpos=wsq*pos
neg=(x<0)*1
xneg=x*neg
wneg=w*neg
wxneg=wx*neg
xsqneg=xsq*neg
wsqneg=wsq*neg
wxsq=w*xsq
xwsq=x*wsq
wcub=w^3
xcub=x^3


modelrd1<-lm(y~pos+xpos+xsqpos+neg+xneg+xsqneg-1)$coefficients
modelrd2<-lm(y~pos+wpos+xpos+xsqpos+neg+wneg+xneg+xsqneg-1)$coefficients
modelrd3<-lm(y~w+x+wx+xsq+wsq+wxsq+xwsq+wcub+xcub)$coefficients
modelfs1<-glm(t~pos+xpos+xsqpos+neg+xneg+xsqneg-1,family=binomial(link="probit"))$coefficients
modelfs2<-glm(t~pos+wpos+xpos+wxpos+xsqpos+wsqpos+neg+wneg+xneg+wxneg+xsqneg+wsqneg-1,family=binomial(link="probit"))$coefficients
save(modelrd1,modelrd2,modelfs1,modelfs2,file="./simulation/DGP_MonoTest.RData")
