\documentclass[notitlepage]{article}
\usepackage{Math598}
\usepackage{listings}
\usepackage{numprint}
\usepackage{enumerate}
\usepackage{bbm}
\usepackage{verbatim}
\usepackage{amsfonts,amsmath}
\usepackage{amssymb}

\def\E{\Expect}

\newcommand{\independent}{\perp\mkern-9.5mu\perp}

\lstloadlanguages{R}

\definecolor{keywordcolor}{rgb}{0,0.6,0.6}
\definecolor{delimcolor}{rgb}{0.461,0.039,0.102}
\definecolor{Rcommentcolor}{rgb}{0.101,0.043,0.432}

\lstdefinestyle{Rsettings}{
  basicstyle=\ttfamily,
  breaklines=true,
  showstringspaces=false,
  keywords={if, else, function, theFunction, tmp}, % Write as many keywords
  otherkeywords={},
  commentstyle=\itshape\color{Rcommentcolor},
  keywordstyle=\color{keywordcolor},
  moredelim=[s][\color{delimcolor}]{"}{"},
}

\lstset{basicstyle=\ttfamily, numbers=none, literate={~} {$\sim$}{2}}

\parindent0in

\begin{document}

<<setup, cache=FALSE, echo=FALSE>>=
library(knitr)
# global chunk options
opts_chunk$set(cache=TRUE, autodep=TRUE)
options(scipen=999)
options(repos=c(CRAN="https://cloud.r-project.org/"))
inline_hook <- function (x) {
  if (is.numeric(x)) {
    # ifelse does a vectorized comparison
    # If integer, print without decimal; otherwise print two places
    res <- ifelse(x == round(x),
      sprintf("%6i", x),
      sprintf("%.2f", x)
    )
    paste(res, collapse = ", ")
  }
}
knit_hooks$set(inline = inline_hook)
@

\begin{center}
{\textsclarge{Causal Contrasts from Experimental Study Data}}
\end{center}

In a randomized experimental study, inference concerning the causal effect of binary variable $Z$ on $Y$ can be made by direct comparison of sample averages.  Suppose that $Z \sim Bernoulli(p)$ for $0<p<1$.  Two estimators of the causal contrast
\[
\delta = \E[Y(1) - Y(0)]
\]
can be proposed:
\begin{equation}\label{eq:tilde}
\widetilde{\delta} = \frac{1}{np} \sum_{i=1}^n Z_i Y_i - \frac{1}{n(1-p)} \sum_{i=1}^n (1-Z_i) Y_i = \frac{1}{n p (1-p)} \sum_{i=1}^n (Z_i - p) Y_i
\end{equation}
and
\begin{equation}\label{eq:hat}
\widehat{\delta} = \frac{1}{n \widehat p} \sum_{i=1}^n Z_i Y_i - \frac{1}{n(1-\widehat p)} \sum_{i=1}^n (1-Z_i) Y_i = \frac{1}{n \widehat p (1- \widehat p)} \sum_{i=1}^n (Z_i - \widehat p) Y_i
\end{equation}
where
\[
\widehat p = \frac{1}{n} \sum_{i=1}^n Z_i.
\]
It follows that \eqref{eq:hat} has lower variance than \eqref{eq:tilde}.

\bigskip

\begin{enumerate}[(1)]

\item The variance of $\widetilde{\delta}$ is computed as follows: we have that
\[
\widetilde{\delta} = \frac{1}{np(1-p)} \sum_{i=1}^n (Z_i-p) Y_i
\]
which has variance
\[
\frac{1}{n p^2 (1-p)^2} \Var[(Z-p) Y].
\]
In the experimental study, $Z \independent X$.
\begin{comment}
\begin{align*}
  \Var[(Z-p)Y] & = \E_{X,Z}[(Z-p)^2 \Var_{Y|X,Z}[Y|X,Z]] + \Var_{X,Z}[(Z-p)\E_{Y|X,Z}[Y|X,Z]] \\[6pt]
   & = \E_{X,Z}[(Z-p)^2] \sigma_Y^2 + \Var_{X,Z}[(Z-p)(\alpha X + \delta Z)]\\[6pt]
   & = p(1-p) \sigma_Y^2 + \E_{X,Z}[(Z-p)^2(\alpha X + \delta Z)^2] - \left\{\E_{X,Z}[(Z-p)(\alpha X + \delta Z)] \right\}^2 \\[6pt]
   & = p(1-p) \sigma_Y^2 + \E_{X,Z}[(Z-p)^2(\alpha^2 X^2 + 2 \alpha \delta X + \delta^2 Z^2)] - (\delta p(1-p))^2 \\[6pt]
   & = p(1-p) \sigma_Y^2 + \alpha^2 m_2 p(1-p) + 2 \alpha \delta m_1 p(1-p)^2 + \delta^2 p (1-p)^2 - \delta^2 p^2 (1-p)^2
\end{align*}
where $m_k = \E[X^k]$ and
\[
\E[Z(Z-p)^2] \E[Z (Z-p)^2] = \E[Z^2 (Z-p)^2] = p (1-p)^2.
\]
Therefore
\[
\Var[\widetilde \delta] = \frac{1}{p(1-p)} \left(\sigma_Y^2 + \alpha^2 m_2 + 2 \alpha \delta m_1 (1-p) + \delta^2 (1-p)^2 \right)
\]
\end{comment}

\begin{align*}
  \Var[(Z-p)Y] & = \E_{Z}[(Z-p)^2 \Var_{Y|Z}[Y|Z]] + \Var_{Z}[(Z-p)\E_{Y|Z}[Y|Z]].
\end{align*}
We have that
\[
\E_{Y|Z}[Y|Z] = \E_X[\E_{Y|X,Z}[Y|X,Z]] = \E_X[\mu(X,Z)] = \mu(Z)
\]
say.  Note that
\[
\E_Z[(Z-p) \mu(Z)] = -p(1-p) \mu(0) + p(1-p) \mu(1) = p(1-p) (\mu(1)-\mu(0)) = p(1-p) \delta.
\]
If $\Var_{Y|X,Z}[Y|X,Z] = \sigma^2(X,Z)$, we also have that
\begin{align*}
\Var_{Y|Z}[Y|Z] & = \E_X[\Var_{Y|X,Z}[Y|X,Z]] + \Var_X[\E_{Y|X,Z}[Y|X,Z]] \\[6pt]
& = \E_X[\sigma^2(X,Z)] + \Var_X[\mu(X,Z)] = \sigma^2(Z)
\end{align*}
say. Therefore
\begin{align*}
  \Var[(Z-p)Y] & = \E_{Z}[(Z-p)^2 \sigma^2(Z)] + \Var_{Z}[(Z-p)\mu(Z)] \\[6pt]
  & = \left\{p^2 (1-p) \sigma^2(0) + p (1-p)^2 \sigma^2(1)  \right\} + \left\{\E_{Z}[(Z-p)^2{\mu(Z)}^2] - \{\E_Z[(Z-p) \mu(Z)]\}^2\right\}\\[6pt]
  & = p(1-p) \left\{ p \sigma^2(0) + (1-p) \sigma^2(1)\right\} \\[6pt]
  & \qquad +  p^2 (1-p) \{\mu(0)\}^2 + p (1-p)^2 \{\mu(0)\}^2 - p^2(1-p)^2 (\mu(1)-\mu(0))^2\\[6pt]
  & = p(1-p) \left\{ p \sigma^2(0) + (1-p) \sigma^2(1 )\right\} + p(1-p) (p \mu(0) + (1-p) \mu(1))^2
\end{align*}
Therefore
\[
\Var[\widetilde \delta] = \frac{1}{n p (1-p)} \left( p \sigma^2(0) + (1-p) \sigma^2(1) + (p \mu(0) + (1-p) \mu(1))^2 \right).
\]

\pagebreak

\item The variance of $\widehat{\delta}$
\[
\widehat{\delta} = \frac{1}{n \widehat p(1-\widehat p)} \sum_{i=1}^n (Z_i-\widehat p) Y_i
\]
can be computed by noting that the estimator arises as the WLS estimator of $\delta$ from the model
\[
Y_i = \beta_0 + Z_i \delta + \varepsilon_i
\]
where
\[
\Var_{\varepsilon|Z}[\varepsilon|Z=z] = \sigma^2(z).
\]
Under homoscedasticity, $\sigma^2(z) \equiv \sigma^2$.

\medskip

For large $n$, the variance covariance matrix of the WLS estimators is known to be the inverse of the expectation of
\[
\frac{n}{\sigma^2(Z)}
\begin{bmatrix}
   1 & Z \\
   Z & Z^2
\end{bmatrix}.
\]
Taking expectations (elementwise) we have
\[
n
\begin{bmatrix}
   \dfrac{(1-p)}{\sigma^2(0)} + \dfrac{p}{\sigma^2(1)} & \dfrac{p}{\sigma^2(1)} \\[12pt]
   \dfrac{p}{\sigma^2(1)} & \dfrac{p}{\sigma^2(1)}
\end{bmatrix}
\]
which has inverse
\[
\frac{1}{n} \frac{\sigma^2(0) \sigma^2(1)}{p(1-p)} 
\begin{bmatrix}
   \dfrac{p}{\sigma^2(1)} & -\dfrac{p}{\sigma^2(1)} \\[12pt]
   -\dfrac{p}{\sigma^2(1)} & \dfrac{(1-p)}{\sigma^2(0)} + \dfrac{p}{\sigma^2(1)} 
\end{bmatrix}
\]
Picking off the element (2,2) element of this matrix, we find that for large $n$
\[
\Var[\widehat \delta] \bumpeq \frac{1}{n p (1-p)} (p\sigma^2(0)+(1-p) \sigma^2(1)).
\]
\end{enumerate}

We therefore have from (1) that
\[
\Var[\widetilde \delta] = \frac{1}{n p (1-p)} \left( p \sigma^2(0) + (1-p) \sigma^2(1) + (p \mu(0) + (1-p) \mu(1))^2 \right) \geq \Var[\widehat \delta].
\]
Most commonly we will assume homoscedasticity where $\sigma^2(z) \equiv \sigma^2$, in which case
\begin{align*}
  \Var[\widetilde \delta] & = \dfrac{1}{n p (1-p)} \left( \sigma^2 + (p \mu(0) + (1-p) \mu(1))^2 \right) \\[6pt]
  \Var[\widehat \delta] & = \dfrac{1}{n p (1-p)} \sigma^2
\end{align*}

\pagebreak

\textbf{Simulation 1:} In the simulation below, we assume that
\[
X \sim Normal(\mu_X,\sigma_X^2)
\]
with $\mu_X = 1$ and $\sigma_X = 1$, and
\[
Y | X = x, Z = z \sim Normal( \alpha x + \delta z, \sigma^2)
\]
with $\sigma^2(x,z) \equiv \sigma^2 = 1$, and
\[
\mu(x,z) = \E_{Y|X,Z}[Y|X=x,Z=z] = \alpha x + \delta z
\]
with $\alpha = 0.5$ and $\delta = 2$, and hence
\[
\E[Y(1) - Y(0)] = \E_{Y|X,Z}[Y|X=x,Z=1] -  \E_{Y|X,Z}[Y|X=x,Z=0] = \delta .
\]
We also have
\[
\mu(z) = \E_X[\mu(X,z)] = \E_X[\alpha X + \delta z] = \alpha \E_X[X]  + \delta z = \alpha \mu_X + \delta z
\]
and that
\begin{align*}
\sigma^2(z) & = \E_X[\sigma^2(X,z)] + \Var_X[\mu(X,z)] \\[6pt]
& = \sigma^2 + \Var_X[\alpha X + \delta z] \\[6pt]
& = \sigma^2 + \alpha^2 \Var_X[X]\\[6pt]
& = \sigma^2 + \alpha^2 \sigma_X^2
\end{align*}

\medskip

We study how the variances change as $n$ increases.  We choose $p=0.2$.
<<exp5,comment='+', fig.width=7.5, fig.height=4.5, fig.align='center'>>=
set.seed(23)
nreps<-10000;nvec<-seq(50,2000,by=50)
p<-0.2;delta<-2;al<-0.5;sig<-1
muX<-1;sigX<-1
mu0<-al*muX; mu1<-mu0+delta
ests.mat<-array(0,c(length(nvec),nreps,2))
for(i in 1:length(nvec)){
    n<-nvec[i]
    for(irep in 1:nreps){
        X<-rnorm(n,muX,sigX)
        Z<-rbinom(n,1,p)
        Y<-rnorm(n,delta*Z+al*X,sig)
        p.hat<-mean(Z)
        ests.mat[i,irep,1]<-sum(Z*Y)/(n*p)-sum((1-Z)*Y)/(n*(1-p))
        ests.mat[i,irep,2]<-sum(Z*Y)/(n*p.hat)-sum((1-Z)*Y)/(n*(1-p.hat))
    }
}
#Variance of delta tilde
true.var.tilde<-((sig^2+al^2*sigX^2+(p*mu0+(1-p)*mu1)^2)/(p*(1-p)))
true.var.tilde

#Variance of delta hat
true.var.hat<-(sig^2+al^2*sigX^2)/(p*(1-p))
true.var.hat
#Efficiency
true.var.tilde/true.var.hat
@
<<exp6,comment='+', fig.width=7.5, fig.height=4.9, fig.align='center',echo=FALSE>>=
par(mar=c(4,2,3,1),pch=19)
boxplot(t(ests.mat[,,1]),ylim=range(0.5,3),cex=0.5,names=nvec,cex.axis=0.75);
abline(h=delta,lty=2,col='red')
title(expression(paste('Variability of ',widetilde(delta))))
boxplot(t(ests.mat[,,2]),ylim=range(0.5,3),cex=0.5,names=nvec,cex.axis=0.75);
abline(h=delta,lty=2,col='red')
title(expression(paste('Variability of ',widehat(delta))))
@

<<exp7,comment='+', fig.width=7.5, fig.height=5, fig.align='center'>>=
v1<-apply(ests.mat[,,1],1,var)
#Variance for tilde delta
nvec*v1
par(mar=c(4,4,3,1))
plot(nvec,nvec*v1,type='l',xlab='n',ylab='Variance')
points(nvec,nvec*v1,pch=19,cex=0.5)
abline(h=true.var.tilde,col='red',lty=2)
title(expression(paste('Variance of ',sqrt(n)*widetilde(delta))))
@

<<exp8,comment='+', fig.width=7.5, fig.height=5, fig.align='center'>>=
v2<-apply(ests.mat[,,2],1,var)
#Variance for delta hat
nvec*v2
par(mar=c(4,4,3,1))
plot(nvec,nvec*v2,type='l',xlab='n',ylab='Variance')
points(nvec,nvec*v2,pch=19,cex=0.5)
abline(h=true.var.hat,col='red',lty=2)
title(expression(paste('Variance of ',sqrt(n)*widehat(delta))))
@

<<exp9,comment='+', fig.width=7.5, fig.height=5.0, fig.align='center',echo=FALSE>>=
eff<-apply(ests.mat[,,1],1,var)/apply(ests.mat[,,2],1,var)
par(mar=c(4,4,2,1))
plot(nvec,eff,type='l',xlab='n',ylab='Efficiency')
points(nvec,eff,pch=19,cex=0.5)
abline(h=true.var.tilde/true.var.hat,col='red',lty=2)
title('Efficiency')
@
The efficiency is the ratio of variances of the two estimators: $\Var[\widetilde \delta]/\Var[\widehat \delta]$.

\pagebreak

\textbf{Simulation 2:} In the simulation below, we assume that
\[
X \sim Normal(\mu_X,\sigma_X^2)
\]
with $\mu_X = 1$ and $\sigma_X = 1$, and
\[
Y | X = x, Z = z \sim Normal( \alpha x + \delta z, \sigma^2 (x,z))
\]
with $\sigma^2(x,z) = \sigma^2(1 + z)$, and
\[
\mu(x,z) = \E_{Y|X,Z}[Y|X=x,Z=z] = \alpha x + \delta z
\]
with $\alpha = 0.5$ and $\delta = 2$, and hence
\[
\E[Y(1) - Y(0)] = \E_{Y|X,Z}[Y|X=x,Z=1] -  \E_{Y|X,Z}[Y|X=x,Z=0] = \delta .
\]
We also have
\[
\mu(z) = \E_X[\mu(X,z)] = \E_X[\alpha X + \delta z] = \alpha \E_X[X]  + \delta z = \alpha \mu_X + \delta z
\]
and that
\begin{align*}
\sigma^2(z) & = \E_X[\sigma^2(X,z)] + \Var_X[\mu(X,z)] \\[6pt]
& = \sigma^2 (1 + z) + \Var_X[\alpha X + \delta z] \\[6pt]
& = \sigma^2 (1 + z) + \alpha^2 \Var_X[X]\\[6pt]
& = \sigma^2 (1 + z) + \alpha^2 \sigma_X^2
\end{align*}

\medskip

We study how the variances change as $n$ increases.  We choose $p=0.2$.
<<exp10,comment='+', fig.width=7.5, fig.height=4.5, fig.align='center'>>=
set.seed(23)
nreps<-10000;nvec<-seq(50,2000,by=50)
p<-0.2;delta<-2;al<-0.5;sig<-1
muX<-1;sigX<-1
mu0<-al*muX; mu1<-mu0+delta
sigs<-c(0,0)
sigs[1]<-sqrt(sig^2*(1)+al^2*sigX^2)
sigs[2]<-sqrt(sig^2*(1+1)+al^2*sigX^2)
ests.mat<-array(0,c(length(nvec),nreps,2))
for(i in 1:length(nvec)){
    n<-nvec[i]
    for(irep in 1:nreps){
        X<-rnorm(n,muX,sigX)
        Z<-rbinom(n,1,p)
        Y<-rnorm(n,delta*Z+al*X,sigs[Z+1])
        p.hat<-mean(Z)
        ests.mat[i,irep,1]<-sum(Z*Y)/(n*p)-sum((1-Z)*Y)/(n*(1-p))
        ests.mat[i,irep,2]<-sum(Z*Y)/(n*p.hat)-sum((1-Z)*Y)/(n*(1-p.hat))
    }
}
#Variance of delta tilde
true.var.tilde<-(((1-p)*sigs[2]^2+p*sigs[1]^2)+al^2*sigX^2+(p*mu0+(1-p)*mu1)^2)/(p*(1-p))
true.var.tilde

#Variance of delta hat
true.var.hat<-(((1-p)*sigs[2]^2+p*sigs[1]^2)+al^2*sigX^2)/(p*(1-p))
true.var.hat
#Efficiency
true.var.tilde/true.var.hat
@
<<exp11,comment='+', fig.width=7.5, fig.height=4.9, fig.align='center',echo=FALSE>>=
par(mar=c(4,2,3,1),pch=19)
boxplot(t(ests.mat[,,1]),ylim=range(0.5,3),cex=0.5,names=nvec,cex.axis=0.75);
abline(h=delta,lty=2,col='red')
title(expression(paste('Variability of ',widetilde(delta))))
boxplot(t(ests.mat[,,2]),ylim=range(0.5,3),cex=0.5,names=nvec,cex.axis=0.75);
abline(h=delta,lty=2,col='red')
title(expression(paste('Variability of ',widehat(delta))))
@

<<exp12,comment='+', fig.width=7.5, fig.height=5, fig.align='center'>>=
v1<-apply(ests.mat[,,1],1,var)
#Variance for tilde delta
nvec*v1
par(mar=c(4,4,3,1))
plot(nvec,nvec*v1,type='l',xlab='n',ylab='Variance')
points(nvec,nvec*v1,pch=19,cex=0.5)
abline(h=true.var.tilde,col='red',lty=2)
title(expression(paste('Variance of ',sqrt(n)*widetilde(delta))))
@

<<exp13,comment='+', fig.width=7.5, fig.height=5, fig.align='center'>>=
v2<-apply(ests.mat[,,2],1,var)
#Variance for delta hat
nvec*v2
par(mar=c(4,4,3,1))
plot(nvec,nvec*v2,type='l',xlab='n',ylab='Variance')
points(nvec,nvec*v2,pch=19,cex=0.5)
abline(h=true.var.hat,col='red',lty=2)
title(expression(paste('Variance of ',sqrt(n)*widehat(delta))))
@

<<exp14,comment='+', fig.width=7.5, fig.height=5.0, fig.align='center',echo=FALSE>>=
eff<-apply(ests.mat[,,1],1,var)/apply(ests.mat[,,2],1,var)
par(mar=c(4,4,2,1))
plot(nvec,eff,type='l',xlab='n',ylab='Efficiency')
points(nvec,eff,pch=19,cex=0.5)
abline(h=true.var.tilde/true.var.hat,col='red',lty=2)
title('Efficiency')
@

\pagebreak

\textbf{Simulation 3:} In the simulation below, we assume that
\[
X \sim Normal(\mu_X,\sigma_X^2)
\]
with $\mu_X = 1$ and $\sigma_X = 1$, and
\[
Y | X = x, Z = z \sim Normal( \alpha x + \delta z, \sigma^2 (x,z))
\]
with $\sigma^2(x,z) = \sigma^2(x^2 + z)$, and
\[
\mu(x,z) = \E_{Y|X,Z}[Y|X=x,Z=z] = \alpha x + \delta z
\]
with $\alpha = 0.5$ and $\delta = 2$, and hence
\[
\E[Y(1) - Y(0)] = \E_{Y|X,Z}[Y|X=x,Z=1] -  \E_{Y|X,Z}[Y|X=x,Z=0] = \delta .
\]
We also have
\[
\mu(z) = \E_X[\mu(X,z)] = \E_X[\alpha X + \delta z] = \alpha \E_X[X]  + \delta z = \alpha \mu_X + \delta z
\]
and that
\begin{align*}
\sigma^2(z) & = \E_X[\sigma^2(X,z)] + \Var_X[\mu(X,z)] \\[6pt]
& = \sigma^2 (\E_X[X^2] + z) + \Var_X[\alpha X + \delta z] \\[6pt]
& = \sigma^2 (\mu_X^2+\sigma_X^2 + z) + \alpha^2 \Var_X[X]\\[6pt]
& = \sigma^2 (\mu_X^2+\sigma_X^2 + z) + \alpha^2 \sigma_X^2
\end{align*}

\medskip

We study how the variances change as $n$ increases.  We choose $p=0.2$.
<<exp15,comment='+', fig.width=7.5, fig.height=4.5, fig.align='center'>>=
set.seed(23)
nreps<-10000;nvec<-seq(50,2000,by=50)
p<-0.2;delta<-2;al<-0.5;sig<-1
muX<-1;sigX<-1
mu0<-al*muX; mu1<-mu0+delta
sigs<-c(0,0)
sigs[1]<-sqrt(sig^2*(muX^2+sigX^2)+al^2*sigX^2)
sigs[2]<-sqrt(sig^2*(muX^2+sigX^2+1)+al^2*sigX^2)
ests.mat<-array(0,c(length(nvec),nreps,2))
for(i in 1:length(nvec)){
    n<-nvec[i]
    for(irep in 1:nreps){
        X<-rnorm(n,muX,sigX)
        Z<-rbinom(n,1,p)
        Y<-rnorm(n,delta*Z+al*X,sigs[Z+1])
        p.hat<-mean(Z)
        ests.mat[i,irep,1]<-sum(Z*Y)/(n*p)-sum((1-Z)*Y)/(n*(1-p))
        ests.mat[i,irep,2]<-sum(Z*Y)/(n*p.hat)-sum((1-Z)*Y)/(n*(1-p.hat))
    }
}
#Variance of delta tilde
true.var.tilde<-(((1-p)*sigs[2]^2+p*sigs[1]^2)+al^2*sigX^2+(p*mu0+(1-p)*mu1)^2)/(p*(1-p))
true.var.tilde

#Variance of delta hat
true.var.hat<-(((1-p)*sigs[2]^2+p*sigs[1]^2)+al^2*sigX^2)/(p*(1-p))
true.var.hat
#Efficiency
true.var.tilde/true.var.hat
@
<<exp16,comment='+', fig.width=7.5, fig.height=4.9, fig.align='center',echo=FALSE>>=
par(mar=c(4,2,3,1),pch=19)
boxplot(t(ests.mat[,,1]),ylim=range(0.5,3),cex=0.5,names=nvec,cex.axis=0.75);
abline(h=delta,lty=2,col='red')
title(expression(paste('Variability of ',widetilde(delta))))
boxplot(t(ests.mat[,,2]),ylim=range(0.5,3),cex=0.5,names=nvec,cex.axis=0.75);
abline(h=delta,lty=2,col='red')
title(expression(paste('Variability of ',widehat(delta))))
@

<<exp17,comment='+', fig.width=7.5, fig.height=5, fig.align='center'>>=
v1<-apply(ests.mat[,,1],1,var)
#Variance for tilde delta
nvec*v1
par(mar=c(4,4,3,1))
plot(nvec,nvec*v1,type='l',xlab='n',ylab='Variance')
points(nvec,nvec*v1,pch=19,cex=0.5)
abline(h=true.var.tilde,col='red',lty=2)
title(expression(paste('Variance of ',sqrt(n)*widetilde(delta))))
@

<<exp18,comment='+', fig.width=7.5, fig.height=5, fig.align='center'>>=
v2<-apply(ests.mat[,,2],1,var)
#Variance for delta hat
nvec*v2
par(mar=c(4,4,3,1))
plot(nvec,nvec*v2,type='l',xlab='n',ylab='Variance')
points(nvec,nvec*v2,pch=19,cex=0.5)
abline(h=true.var.hat,col='red',lty=2)
title(expression(paste('Variance of ',sqrt(n)*widehat(delta))))
@

<<exp19,comment='+', fig.width=7.5, fig.height=5.0, fig.align='center',echo=FALSE>>=
eff<-apply(ests.mat[,,1],1,var)/apply(ests.mat[,,2],1,var)
par(mar=c(4,4,2,1))
plot(nvec,eff,type='l',xlab='n',ylab='Efficiency')
points(nvec,eff,pch=19,cex=0.5)
abline(h=true.var.tilde/true.var.hat,col='red',lty=2)
title('Efficiency')
@
\end{document}