rm(list=ls())

library(pBrackets)

## simulated example

set.seed(776)
pdf("RDiscon_sim.pdf")
n = 50
X = runif(n,0,4)
Y = (1/3)*X + rnorm(n,0,.2)  
Y[X>2] = Y[X>2] + 1        
plot(X,Y,pch=21,bg="lightpink",cex=2,xlim=c(0,4),ylim=c(-.5,3),xlab="",ylab="",
                xaxt="n",yaxt="n")
segments(0,0,2,2/3,col="purple",lwd=3.5)
segments(2,2/3+1,4,4/3+1,col="purple",lwd=3.5)
abline(v=2,col="red",lwd=3)
brackets(2,2/3,2,2/3+1,curvature=1,lwd=3)
title(ylab="Y", line=1, cex.lab=2)
axis(side=1, at=2,"c",tck=0, col.axis="red",cex.axis=2,cex.lab=2)
mtext("x",side=1,line=1,cex=2,adj=.9)

dev.off()




#### data analysis

kernel = function(x,y,grid,h){
     m = length(grid)
     out = rep(0,m)
     for(i in 1:m){
          w = dnorm(grid[i],x,h)
          out[i] = sum(y*w)/sum(w)
          }

     return(out)
     }


Local_Linear = function(x,y,grid,h){
     m = length(grid)
     out = rep(0,m)
     for(i in 1:m){
          w = dnorm(grid[i],x,h)
          tmp = lm(y ~ x,weights=w)
          out[i] = tmp$coef[1] + grid[i]*tmp$coef[2]
          }

     return(out)
     }


pdf("RDiscon_example.pdf")

library(rddtools)
data(house)
x = house$x
y = house$y

h = .3

##left
grid0 = seq(-1,0,length=100)
xx = x[x<0]
yy = y[x<0]
mu0 = Local_Linear(xx,yy,grid0,h)
w = dnorm(grid0[100],x,h)
tmp = lm(y ~ x,weights=w)
beta00 = tmp$coef[1]
se0 = summary(tmp)$coefficients[1,2]


##rigth 
grid1 = seq(0,1,length=100)
xx = x[x>0]
yy = y[x>0]
mu1 = Local_Linear(xx,yy,grid1,h)
w = dnorm(grid1[1],x,h)
tmp = lm(y ~ x,weights=w)
beta01 = tmp$coef[1]
se1 = summary(tmp)$coefficients[1,2]


##point estimate

print(beta00)
print(se0)
print(beta01)
print(se1)
psi = beta01 - beta00
se = sqrt(se0^2 + se1^2)
CI = c(psi-2*se,psi+2*se)
print(CI)
print(se)


## plot
mu = c(mu0,mu1)
grid = c(grid0,grid1)
plot(grid,mu,type="n",xlab="",ylab="")
lines(grid0,mu0,lwd=4,col="purple")
lines(grid1,mu1,lwd=4,col="purple")
abline(v=0,col="red",lty=4,lwd=3)
points(x,y,pch=".",col="blue")
title(ylab="Current Election", line=1.8, cex.lab=2)
title(xlab="Previous Election", cex.lab=2)

dev.off()