# Install and load required packages
install.packages("ivDiag")
install.packages("haven")
install.packages("AER")
install.packages("sandwich")
install.packages("lmtest")
library(ivDiag)
library(haven)
library(AER)       # For instrumental variable regression
library(sandwich)  # For robust standard errors
library(lmtest)    # For hypothesis tests and robust SEs

# Read the Fulton Fish Market data
fulton <- read_dta("https://github.com/scunning1975/mixtape/raw/master/Fulton.dta")

# Label variables for clarity
# q: Log quantity of whiting sold in pounds
# p: Log average daily price per pound
# Stormy: Instrument for p

# Convert day-of-week variables to factors
fulton$Mon <- as.factor(fulton$Mon)
fulton$Tue <- as.factor(fulton$Tue)
fulton$Wed <- as.factor(fulton$Wed)
fulton$Thu <- as.factor(fulton$Thu)

# OLS regression with robust standard errors
ols_model <- lm(q ~ p + Mon + Tue + Wed + Thu, data = fulton)
ols_robust <- coeftest(ols_model, vcov = vcovHC(ols_model, type = "HC1"))
cat("OLS Regression Results with Robust Standard Errors:\n")
print(ols_robust)

# 2SLS regression using 'Stormy' as an instrument for 'p'
iv_formula <- as.formula("q ~ p + Mon + Tue + Wed + Thu | Stormy + Mon + Tue + Wed + Thu")
iv_model <- ivreg(iv_formula, data = fulton)
iv_robust <- coeftest(iv_model, vcov = vcovHC(iv_model, type = "HC1"))
cat("\n2SLS Regression Results with Robust Standard Errors:\n")
print(iv_robust)

# First-stage regression: p ~ Stormy + controls
first_stage <- lm(p ~ Stormy + Mon + Tue + Wed + Thu, data = fulton)
first_stage_robust <- coeftest(first_stage, vcov = vcovHC(first_stage, type = "HC1"))
cat("\nFirst-Stage Regression Results:\n")
print(first_stage_robust)

# Calculate Olea-Pflueger effective F statistic
effF <- eff_F(data = fulton,
              Y = "q",      # Dependent variable
              D = "p",      # Endogenous variable
              Z = "Stormy", # Instrument
              cl = NULL,    # No clustering
              weights = NULL)  # No weights
cat("\nOlea-Pflueger Effective F Statistic:\n")
print(effF)

# Perform Anderson-Rubin test with confidence intervals
ar_results <- AR_test(data = fulton,
                      Y = "q",
                      D = "p",
                      Z = "Stormy",
                      controls = c("Mon", "Tue", "Wed", "Thu"),
                      CI = TRUE,  # Confidence interval
                      alpha = 0.05)  # 5% significance level
cat("\nAnderson-Rubin Test Results and Confidence Intervals:\n")
print(ar_results)

# Perform complete diagnostic analysis with ivDiag
iv_results <- ivDiag(data = fulton,
                     Y = "q",
                     D = "p",
                     Z = "Stormy",
                     controls = c("Mon", "Tue", "Wed", "Thu"),
                     bootstrap = TRUE,
                     run.AR = TRUE)

# Print and plot the results
cat("\nIV Diagnostics Results:\n")
print(iv_results)
plot_coef(iv_results)