A basic implementation of Bayesian (linear) instrumental variables regression with a single endogenous variable and a single IV. This is just a byproduct of me trying to learn more about Turing.jl and Bayesian estimation, but I'm hoping it's useful to someone as an example/starting point.

BayesianIV.jl

using Plots
using StatsPlots, Turing 
using Distributions, LinearAlgebra, Random, FillArrays

Random.seed!(12334455)

# -------------------------------------------------------
# Models/functions for OLS and IV 
# Only implemented with a single RHS variable -- easy 
# to extend to multiple 
# -------------------------------------------------------
@model function bayesian_linreg(data)
    N = length(data[:,1])
    X = data[:,2]
    
    # Priors for alpha and beta
    alpha ~ Normal(0, 0.2)
    beta ~ Normal(0, 0.5)
    
    # Prior for sigma
    sigma ~ Exponential(1)
    
    # Calculate mu
    mu = alpha .+ X * beta
    
    # Likelihood
    for i in 1:N
        data[i,1] ~ Normal(mu[i], sigma)
    end
end

@model function bayesian_iv_onevar(data, IV, endog_cols)
    Y = data[:,1]
    endog = data[:,endog_cols]

    N = length(Y)

    # Priors for alphaY and alphaX
    alphaY ~ Normal(0, 0.2)
    alphaX ~ Normal(0, 0.2)
    
    # Prior for beta
    beta ~ Normal(0, 0.5)
    
    # Prior for gamma
    gamma ~ Normal(0, 0.5)
    
    # Prior for correlation matrix Rho
    Rho ~ LKJ(2, 1)
    
    # Prior for scale parameters Sigma
    Sigma ~ Exponential(1)
    
    # Calculate muY and muX
    muX = alphaX .+ gamma .* IV
    muY = alphaY .+ beta .* muX
    
    # Calculate covariance matrix
    cov_matrix = LinearAlgebra.Diagonal(FillArrays.Fill(Sigma, 2)) .* Rho .* LinearAlgebra.Diagonal(FillArrays.Fill(Sigma, 2));
    
    # Likelihood for Y and X
    # outcomes = hcat(Y, endog)
    for i in 1:N
        data[i,:] ~ MvNormal([muY[i], muX[i]], cov_matrix)
    end
end

# Simulate a linear IV model
N = 2000
error = randn(N)
IV = randn(N) 
X = 1.3 .*IV + 0.7 .* error;
Y = X .* 0.5 + 0.5 .* error;

# Estimation 
chain_ols = sample(bayesian_linreg([Y X]), Turing.NUTS(), 10_000); # linear regression
chain = sample(bayesian_iv_onevar([Y X], IV, 2), Turing.NUTS(), 10_000); # IV

# Plot chains 
density(chain["beta"], label = "Posterior, IV") # pretty good! we beat endogeneity
density!(chain_ols["beta"], label = "Posterior, OLS") # bad -- endogeneity biases us
vline!([0.5], label = "True value")