How to calculate the variance of an estimator with simulation in R So I have a question in the exam of Statistics that I dont know how to do it:
There is $X\sim \text{Poisson}(\theta)$.
I want to estimate $\exp[-\theta] = P(X=0)$.
I know by the maximum likelihood method that a consistent estimator is $\exp[-\bar{X}]$. 
Now is asked me to use R commands to calculate an approximation of $\text{Var}(\exp[-\bar{X}])$, with $n=20$.
I am really not getting how I should do this if I don't have any knowledge about the parameter $\theta$. How I should simulate the data? 
 A: Here is a short R code comparing the estimators
$$\frac{1}{n}\sum_{t=1}^n \mathbb{I}_0(X_t)$$
(in yellow)
and $$\exp\left\{-\sum_{t=1}^nX_t\big/n\right\}$$
(in brown) for $n=20$ for a range of values of $\theta$
compest<-function(T=1e3,tmin=.1,tmax=10){

 varz=matrix(0,50,2)
 theta=seq(tmin,tmax,le=50)
 for (t in 1:50){
  sampl=matrix(rpois(n=20*T,lambda=theta[t]),ncol=20)
  varz[t,]=c(var(apply(sampl==0,1,mean)),
             var(exp(-apply(sampl,1,mean))))}
 return(varz)}

This shows that the estimator based on $\bar{X}_n$ is leading to a smaller variance than the one based on the frequency of zero draws. 
The continuous curves are the theoretical values of the variances, namely $e^{-\theta}(1-e^{-\theta})/n$ for the Binomial proportion of zero draws and $e^{-2\theta}\theta/n$ for the exponential of the average. (This variance is a delta-method approximation of the exact variance, but the fit is very good!)

Actually, as pointed out by George Henry on my blog, the derivation of the mean and variance of $$\exp\left\{-\sum_{t=1}^nX_t\big/n\right\}$$ is quite manageable: since $n\bar{X}_n$ is a Poisson $\mathscr{P}(n\theta)$ variable
\begin{align*}
\mathbb{E}[\exp\{-\bar{X}_n\}]&=\sum_{i=0}^\infty \exp\{-i/n\}\frac{(n\theta)^i}{i!}\exp\{-n\theta\}\\
&=\sum_{i=0}^\infty \left(\exp\{-1/n\} n\theta \right)^i 
\frac{\exp\{-\theta\}}{i!}\\
&=\exp\left\{-n\theta+n\theta\exp\{-1/n\} \right\}\\
&=\exp\left\{-n\theta[1-\exp\{-1/n\}]\right\}
\end{align*}
and
\begin{align*}
\mathbb{E}[\exp\{-\bar{X}_n\}^2]&=\sum_{i=0}^\infty \exp\{-2i/n\}\frac{(n\theta)^i}{i!}\exp\{-n\theta\}\\
&=\sum_{i=0}^\infty \left(\exp\{-2/n\} n\theta \right)^i 
\frac{\exp\{-\theta\}}{i!}\\
&=\exp\left\{-n\theta+n\theta\exp\{-2/n\} \right\}\\
&=\exp\left\{-n\theta[1-\exp\{-2/n\}]\right\}
\end{align*}
Hence
$$\text{var}(\exp\{-\bar{X}_n\})=\exp\left\{-n\theta[1-\exp\{-2/n\}]\right\}-\exp\left\{-2n\theta[1-\exp\{-1/n\}]\right\}$$
As shown on the plot below, the difference with the approximation is hard to spot!

