A simulation of the method of Instrumental Variables

simIV.R

library(DiagrammeR)
library(DiagrammeRsvg)
library(rsvg)
library(dplyr)
library(dqrng)
library(ggplot2)

g <- grViz("digraph{
                     graph[rankdir = TB]
                     node[shape = none]
                     U[shape=circle]
           { rank = same; W; X; Y }
                     W -> X;
                     U -> X;
                     U -> Y;
                     X -> Y;
                     }")
print(g)
rsvg_png(charToRaw(export_svg(g)), 
         "graph.png", width = 700, height = 700)


simIV <- function(N, alpha_xw){
  #alpha_xw <- 0.2
  alpha_xu <- 0.5
  alpha_yx <- 0.3
  alpha_yu <- 0.3
  #standadize to var=1
  sig_x = sqrt(1-(alpha_xw^2+alpha_xu^2))
  sig_y = sqrt(1-(alpha_yx^2+alpha_yu^2))
  W <- dqrnorm(N)
  U <- dqrnorm(N)
  X <- alpha_xw*W  + alpha_xu*U + dqrnorm(N, sd = sig_x)
  Y <- alpha_yx*X + alpha_yu*U + dqrnorm(N, sd = sig_y)
  betahat = as.vector(solve(t(W)%*%X, t(W)%*%Y))
  return(betahat)
}


rep_simIV <- function(iter, N, alpha_xw){
  replicate(iter, simIV(N = N, alpha_xw = alpha_xw)) 
}


condf = expand.grid(N=c(100, 500, 1000),
                    alpha_xw = seq(-0.6,0.6,by=0.05))
dim(condf)

system.time({
  set.seed(20251013)
  resdf = rowwise(condf) %>% 
    reframe(N=N, alpha_xw = alpha_xw,
            estimates = rep_simIV(10000, N, alpha_xw))  
})
#   user  system elapsed 
# 17.889   0.175  18.227 
p = ggplot(resdf, aes(x=factor(N), y = estimates))+
  geom_violin()+
  geom_hline(yintercept = 0.3, linetype=2) + 
  stat_summary(fun.data = mean_sdl, fun.args = list(mult=1))+
  facet_wrap( ~ round(alpha_xw,2), scales="free_y")+
  labs(x="N") + theme_bw(12)
print(p)
ggsave(filename = "violin.png", plot = p, width = 10, height = 10)