```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Part 1: Simulations on Cauchy

Suppose $\boldsymbol{X}=\left(X_{1}, \ldots, X_{n}\right)$ is an i.i.d. sample from the shifted Cauchy distribution with density
$$
f(x \mid \theta)=\frac{1}{\pi\left(1+(x-\theta)^{2}\right)}, \quad x \in \mathbb{R}
$$

Our goal is to compare the following 4 estimators of the parameter $\theta$. 

- Sample mean
$$
\hat{\theta}_{n}^{(1)}=\bar{X}_{n}=\frac{1}{n} \sum_{i=1}^{n} X_{i}
$$

- Sample median
$$
\hat{\theta}_{\mathrm{n}}^{(2)}=M_{\mathrm{n}}=\frac{1}{2}\left(X_{(k)}+X_{(k+1)}\right)
$$

- Modified sample mean
$$
\hat{\theta}_{n}^{(3)}=M_{n}+\left.\frac{2}{n} \cdot \frac{\partial \ell}{\partial \theta}\right|_{\theta=M_{n}}
$$
where $\ell$ is the log-likelihood function. 

- Maximum likelihood estimator (MLE) $\hat{\theta}_{\mathrm{n}}^{(4)}$ defined by
$$
\ell\left(\hat{\theta}_{n}^{(4)} \mid X\right)=\max _{\theta \in \mathbb{R}} \ell(\theta \mid \boldsymbol{X})
$$
where $\ell$ is the log-likelihood function. 

1. Derive the likelihood function and log-likelihood function. 
2. Simulate data from Cauchy distribution with location 5, and scale 1. 
3. Choose your number of simulations. 
4. Verify consistency of the estimators. There are different approaches. You can samples the data sequentially and plot the sequence of the results as a function of $n$.  What do you observe if it is a consistent estimator? The second approach is to use only the representative increasing values of the sample size. e.g. use $n = 10, 50, 100, 200, \dots, 1000$ and what do you observe?
5. Calculate the mean square error of the estimators.
6. Calculate the coverage probability of the estimators. Calculate $\mathbf{P}_{\theta}\left(\left|\hat{\theta}_{n}-\theta\right| \leq \varepsilon\right)$, for $\varepsilon=0.1$, and $\theta=5$.


<!-- # Part 2: Applying bootstrap  -->

<!-- The following code generates $\left(X_{i}, Y_{i}\right)$ pairs. -->

<!-- ```{r} -->
<!-- library(MASS) -->
<!-- generate_pairs <- function(n) { -->
<!--   # Generate n pairs of financial returns. -->
<!--   muX <- 2 -->
<!--   muY <- -1 -->
<!--   CovMx <- matrix(c(1, -.25, -.25, 2), nrow = 2) -->
<!--   data <- mvrnorm(n = 100, mu = c(muX,muY), Sigma = CovMx) -->
<!--   return(data.frame('X' = data[, 1],  -->
<!--                     'Y' = data[, 2])) -->
<!-- } -->
<!-- ``` -->

<!-- ```{r} -->
<!-- fin_pairs <- generate_pairs( 100 ); # Generate 100 (X,Y) pairs. -->
<!-- head(fin_pairs) -->
<!-- ``` -->

<!-- We are interested in  -->
<!-- $$ -->
<!-- \hat{\alpha}=\frac{\hat{\sigma}_{Y}^{2}-\hat{\sigma}_{X Y}}{\hat{\sigma}_{X}^{2}+\hat{\sigma}_{Y}^{2}-2 \hat{\sigma}_{X Y}}  -->
<!-- $$ -->

<!-- ```{r} -->
<!-- Sigmahat <- cov(fin_pairs)  -->
<!-- Sigmahat -->
<!-- sigma2hatXX <- Sigmahat[1,1] -->
<!-- sigma2hatYY <- Sigmahat[2,2] -->
<!-- sigmahatXY <- Sigmahat[1,2] -->
<!-- ``` -->

<!-- The $\hat{\alpha}$ is -->

<!-- ```{r} -->
<!-- alphahat <- (sigma2hatYY - sigmahatXY)/(sigma2hatXX + sigma2hatYY -2*sigmahatXY) -->
<!-- alphahat -->
<!-- ``` -->

<!-- While the true value of alpha is -->

<!-- ```{r} -->
<!-- sigma2XX <- 1 -->
<!-- sigma2YY <- 2 -->
<!-- sigmaXY <- -0.25 -->
<!-- alpha_true <-(sigma2YY - sigmaXY)/(sigma2XX + sigma2YY -2*sigmaXY) -->
<!-- alpha_true -->
<!-- ``` -->


<!-- Now, again, we're going to resample with replacement from our data, and compute our statistic $\hat{\alpha}$ on each resample. The hope is that these resampled versions of the statistic will resemble the distribution of the statistic evaluated on the original data. -->

<!-- 1. Create a function to compute alphahat from a given data set. -->

<!-- 2. Resample the data $B=200$ times, evaluating $\hat{\alpha}$ on each resample. Then, we'll use those resampled values to estimate the variance. -->

<!-- 3. Create the confidence interval at the estimate.  -->

<!-- ```{r, echo = FALSE} -->
<!-- # https://bookdown.org/larget_jacob/data-modeling-methods/bootstrap.html -->
<!-- ``` -->