## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## -----------------------------------------------------------------------------
strong_cov_mat <- matrix(
	c(
		 #Y      #X1      #X2      #D1      #D2      #Z1      #Z2      #Z3      #Z4
		 3.8635,  0.0387,  0.0922,  1.4810,  0.4778,  0.2024,  0.1085,  0.0850, -0.0368,	#Y
		 0.0387,  0.0832, -0.0029,  0.0415, -0.0014,  0.0025,  0.0009, -0.0004, -0.0058,	#X1
		 0.0922, -0.0029,  0.0848,  0.0068,  0.0462,  0.0029, -0.0009,  0.0028,  0.0027,	#X2
		 1.4810,  0.0415,  0.0068,  1.4345,  0.1949,  0.2184,  0.0562,  0.0464, -0.0155,	#D1
		 0.4778, -0.0014,  0.0462,  0.1949,  0.4042, -0.0034,  0.0909,  0.0875, -0.0079,	#D2
		 0.2024,  0.0025,  0.0029,  0.2184, -0.0034,  0.0851,  0.0005, -0.0043,  0.0000,	#Z1
		 0.1085,  0.0009, -0.0009,  0.0562,  0.0909,  0.0005,  0.0840,  0.0074, -0.0096,	#Z2
		 0.0850, -0.0004,  0.0028,  0.0464,  0.0875, -0.0043,  0.0074,  0.0831, -0.0005,	#Z3
		-0.0368, -0.0058,  0.0027, -0.0155, -0.0079,  0.0000, -0.0096, -0.0005,  0.0827		#Z4
		),
                 ncol = 9)

## -----------------------------------------------------------------------------
weak_cov_mat <- matrix(
	c(
		# Y        X1       X2       D1       D2       Z1       Z2       Z3       Z4
		  0.0795,  0.0009, -0.0003, -0.0022,  0.0053,  0.0003, -0.0009,  0.0002,  0.0010,	#Y
		  0.0009,  0.0832, -0.0029,  0.0020, -0.0036,  0.0025,  0.0009, -0.0004, -0.0058,	#X1
		 -0.0003, -0.0029,  0.0848,  0.0035,  0.0040,  0.0029, -0.0009,  0.0028,  0.0027,	#X2
		 -0.0022,  0.0020,  0.0035,  0.0866, -0.0027,  0.0001,  0.0030, -0.0013, -0.0041,	#D1
		  0.0053, -0.0036,  0.0040, -0.0027,  0.0840, -0.0037,  0.0020, -0.0002, -0.0068,	#D2
		  0.0003,  0.0025,  0.0029,  0.0001, -0.0037,  0.0851,  0.0005, -0.0043,  0.0000,	#Z1
		 -0.0009,  0.0009, -0.0009,  0.0030,  0.0020,  0.0005,  0.0840,  0.0074, -0.0096,	#Z2
		  0.0002, -0.0004,  0.0028, -0.0013, -0.0002, -0.0043,  0.0074,  0.0831, -0.0005,	#Z3
		  0.0010, -0.0058,  0.0027, -0.0041, -0.0068,  0.0000, -0.0096, -0.0005,  0.0827	#Z4
		),
                 ncol = 9)


## -----------------------------------------------------------------------------
library(MASS)
set.seed(42)

mu <- rep(0,9)
stddev <- rep(1,9)
strong_instrument_data <- as.data.frame(MASS::mvrnorm(n = 250, mu = mu, Sigma = strong_cov_mat, empirical = FALSE))
names(strong_instrument_data) <- c("Y", "X1", "X2","D1","D2","Z1","Z2","Z3","Z4")
weak_instrument_data <- as.data.frame(MASS::mvrnorm(n = 250, mu = mu, Sigma = weak_cov_mat, empirical = FALSE))
names(weak_instrument_data) <- c("Y", "X1", "X2","D1","D2","Z1","Z2","Z3","Z4")


## -----------------------------------------------------------------------------
library(cragg)
cragg_donald(X=~X1 + X2, # Control Variables
			 D=~D1 + D2, # Treatments
			 Z=~Z1+Z2+Z3+Z4,# Instruments
			 data = strong_instrument_data)

## -----------------------------------------------------------------------------
stock_yogo_test(X=~X1+X2 , # Control Variables
			 D=~D1 + D2 , # Treatments
			 Z=~Z1+Z2+Z3+Z4 ,# Instruments
			 size_bias="bias", #Default
			 B=.05, #Default
			 data = strong_instrument_data)

## -----------------------------------------------------------------------------
stock_yogo_test(X=~X1+X2 , # Control Variables
			 D=~D1 + D2 , # Treatments
			 Z=~Z1+Z2+Z3+Z4 ,# Instruments
			 size_bias="bias", #Default
			 B=.05, #Default
			 data = weak_instrument_data)