\documentclass[11pt]{article}

\author{Jeremy Morris}
\title{Global Warming \\ \small Project 2}

\usepackage{amssymb,amsmath,amsfonts}
\usepackage[dvips]{graphicx}
\addtolength{\hoffset}{-1.05cm}
\addtolength{\textwidth}{2.5cm}
\addtolength{\textheight}{2.5cm}

\newcommand{\bs}[1]{
\boldsymbol{#1}
}

\begin{document}
\maketitle

\section{A confidence set for the mean vector}
\subsection{$t$ and $T^{2}$ Intervals}
Johnson and Wichern outline a confidence set for a specific linear combination of the mean vector $\bar{\bs{x}}$, which is defined as follows.  \cite{text2}
\begin{equation}
\left( \bs{a}'\bar{\bs{x}} - t_{n-1}(\alpha/2)\sqrt{\frac{\bs{a}'\bs{Sa}}{n}}
\;,\quad
\bs{a}'\bar{\bs{x}} + t_{n-1}(\alpha/2)\sqrt{\frac{\bs{a}'\bs{Sa}}{n}} \; \right)
\label{tinterval}
\end{equation}
Where $t_{n-1}(\alpha/2)$ is the upper $100(\alpha/2)^{th}$ percentile of a $t$-distribution with $n-1$ degrees of freedom.  This confidence set holds for a specific choice of $\bs{a}$.  The following confidence set can be used for $(1-\alpha)$ coverage for all choices of $\bs{a}$.
\begin{equation}
\left(
\bs{a}'\bar{\bs{x}} - \sqrt{\frac{p(n-1)}{n(n-p)}F_{p,n-p}(\alpha)\bs{a}'\bs{Sa}}
\;,\quad
\bs{a}'\bar{\bs{x}} + \sqrt{\frac{p(n-1)}{n(n-p)}F_{p,n-p}(\alpha)\bs{a}'\bs{Sa}}
\;\right) \label{T2interval}
\end{equation}
Where $F_{p,n-p}(\alpha)$ is the upper $100(\alpha)$th percentile of the $F$ distribution with $p$ and $n-p$ degrees of freedom.  This confidence interval is based on the $T^{2}$ statistic and should be used when $(1-\alpha)$ joint coverage is desired for several choices of $\bs{a}$.  The consequence of having joint coverage is that the intervals will be larger than in (\ref{tinterval}).

\subsection{The Bonferroni Method}
For a small set of linear combinations of the mean vector, we can use the Bonferroni method to develop more precise confidence intervals.  \cite{text2}  The Bonferroni method is based on the following probabilistic argument.  If we let $C_{i}$ denote the confidence interval for $\bs{a}'\bs{\mu}$ with $P[C_{i} \textrm{ true} ] = 1 - \alpha_{i}$, where $i = 1,2,\cdots,m$.
\begin{eqnarray}
P[\textrm{all } C_{i} \textrm{ true}] & = & 1 - P[\textrm{at least one } C_{i} \textrm{ false} ] \nonumber \\
& \geq & 1 - \sum_{i=1}^{m}P[C_{i} \textrm{ false}] = 1 - \sum_{i=1}^{m}(1 - P[C_{i} \textrm{ true}]) \nonumber \\
& = & 1 - (\alpha_{1} + \alpha_{2} + \cdots + \alpha_{m})
\end{eqnarray}
This means that we can use the more precise confidence intervals in (\ref{tinterval}) and vary $\alpha$ based on the confidence set we desire.  For example, if we would like to have a $0.95$ confidence set for all means in the temperature data of this project, we would use the linear combinations $\bs{a}_{i}'\bar{\bs{x}}$, with $\bs{a}_{i}' = [0,0,\ldots,a_{i},\ldots,0]'$, where $a_{i} = 1$, and let $\alpha_{i} = \alpha/12$ so that we have the following intervals.
\begin{equation}
\begin{array}{ccccc}
\bar{x}_{1} - t_{214}(\alpha/24)\sqrt{s_{1,1}/215} & \leq & \mu_{1} & \leq & \bar{x}_{1} + t_{214}(\alpha/24)\sqrt{s_{1,1}/215} \\
\bar{x}_{2} - t_{214}(\alpha/24)\sqrt{s_{2,2}/215} & \leq & \mu_{1} & \leq & \bar{x}_{2} + t_{214}(\alpha/24)\sqrt{s_{2,2}/215} \\
&&\vdots&& \\
\bar{x}_{12} - t_{214}(\alpha/24)\sqrt{s_{12,12}/215} & \leq & \mu_{12} & \leq & \bar{x}_{12} + t_{214}(\alpha/24)\sqrt{s_{12,12}/215} \\
\end{array}
\end{equation}

\subsection{Comparison of the $T^{2}$ and Bonferroni intervals}
Notice that the only difference between the intervals in (\ref{tinterval}) and (\ref{T2interval}) is the multiplier corresponding to the critical value.  For (\ref{tinterval}), we use the $t$ distribution and the $F$ distribution for (\ref{T2interval}).  Table (\ref{compare}) shows the multipliers for all three types of intervals used.  
\begin{table}[ht]
\begin{center}
\begin{tabular}{|l|l|l|} \hline
Type & Source & Value \\ \hline
$t$ interval & $t_{214}(0.05/2)$ & 1.971 \\ \hline
Bonferroni $t$ interval & $t_{214}(0.05/24)$ & 2.896 \\ \hline
$T^{2}$ interval & $\sqrt{\frac{p(n-1)}{n-p}F_{12,203}(0.05)}$ & 4.772 \\ \hline
\end{tabular}
\end{center}
\caption{Critical Values} \label{compare}
\end{table}
Notice that the value for the $t$ interval is the smallest, this will give the most precise interval if only one confidence interval is desired.  However, this confidence interval will not give $(1-\alpha)$ confidence for multiple intervals.  Notice also that the multiplier for the $T^{2}$ interval is roughly 60\% larger than the multiplier using the $t$ interval with the Bonferroni method.  Any interval using the $T^{2}$ interval will be twice the length of an interval using the Bonferroni method.  The conclusion reached here is that the most precise confidence intervals will use the $t$-statistic along with the Bonferroni method.

\section{Changing Mean}
In order to test if the mean vector is staying the same with in the data set, we will derive the generalized likelihood ratio test for the following.
\begin{equation}
H_{0} : \bs{\mu}_{1} = \bs{\mu}_{2} = \cdots = \bs{\mu}_{k} \neq \bs{\mu}_{k+1} = \cdots = \bs{\mu}_{n} \qquad H_{a} : \bs{\mu}_{1} = \bs{\mu}_{2} = \cdots = \bs{\mu}_{n} = \bs{\mu}
\label{increase}
\end{equation}
This will test for a single change point within the mean vector against the alternative that there is no change.  In order to use the generalized likelihood ratio test, we will derive the likelihood functions for each hypothesis and determine the MLE's for the parameters.  If we assume that our data comes from a multivariate normal distribution, then our data $\bs{X}_{1}, \bs{X}_{2}, \ldots, \bs{X}_{n}$ has the following density function.
\begin{eqnarray}
f(\bs{x};\bs{\mu},\bs{\Sigma}) & = & \frac{1}{(2\pi)^{p/2}|\bs{\Sigma}|^{1/2}} e^{-\frac{1}{2}(\bs{x} - \bs{\mu})'\bs{\Sigma}^{-1}(\bs{x} - \bs{\mu})}
\end{eqnarray}
The derivation of the likelihood function for the null hypothesis follows.
\begin{eqnarray}
L(\bs{\mu}_{1},\bs{\mu}_{2},k,\bs{\Sigma}) & = & \prod_{i=1}^{k} f(\bs{x}_{i};\bs{\mu}_{1},\bs{\Sigma}) \prod_{j=k+1}^{n} f(\bs{x}_{j};\bs{\mu}_{2},\bs{\Sigma}) \label{eq.null} \\
 & = & L_{1}(\bs{\mu}_{1},\bs{\Sigma})\;L_{2}(\bs{\mu}_{2},\bs{\Sigma}) \label{eq.null.llh}
\end{eqnarray}
Where in (\ref{eq.null.llh}), each likelihood function is defined so that $L_{1}(\bs{\mu}_{1},\bs{\Sigma})$ is taken over the first $k$ observations and $L_{2}(\bs{\mu}_{2},\bs{\Sigma})$ is taken over the last $n - k$ observations.  Then this problem reduces to exercise (6.11) from our text \cite{text2}.  This means that we have the following maximum likelihood estimates.
\begin{eqnarray}
\hat{\bs{\mu}}_{1} & = & \bar{\bs{x}}_{1}\\
\hat{\bs{\mu}}_{2} & = &  \bar{\bs{x}}_{2} \\
\hat{\bs{\Sigma}} & = &  \frac{1}{n}\left[ (k - 1)\bs{S}_{1} + (n - k - 1)\bs{S}_{2} \right]
\end{eqnarray}
Where
\begin{eqnarray}
\bs{S_{1}} & = & \sum_{i=1}^{k}(\bs{x}_{i} - \bar{\bs{x}}_{1})(\bs{x}_{i} - \bar{\bs{x}}_{1})' \\
\bs{S}_{2} & = & \sum_{j=k+1}^{n}(\bs{x}_{j} - \bar{\bs{x}}_{2})(\bs{x}_{j} - \bar{\bs{x}}_{2})'
\end{eqnarray}
There is no closed form maximum likelihood estimate for the parameter $k$.  Since we are looking for the maximum of the likelihood function, we can simply run through all choices of $k$ in the data set for the function $L(\hat{\bs{\mu}}_{1},\hat{\bs{\mu}}_{2},k,\hat{\bs{\Sigma}})$ to determine where the maximum is reached, the particular value of $k$ will be our MLE $\hat{k}$.

The likelihood function is more easily defined for the alternative hypothesis.  Instead of (\ref{eq.null}), we have the following.
\begin{eqnarray}
L(\bs{\mu},\bs{\Sigma}) & = & \prod_{i=1}^{n} f(\bs{X}_{i};\bs{\mu},\bs{\Sigma}) \label{alt}
\end{eqnarray}
With the the maximum reached at the MLE's that follow.
\begin{eqnarray}
\hat{\bs{\mu}} & = &  \bar{\bs{x}} \label{alt.mu} \\
\hat{\bs{\Sigma}} & = &  \left(\frac{n-1}{n}\right)\bs{S} \label{alt.Sigma}
\end{eqnarray}
and
\begin{eqnarray}
\bs{S} & = & \sum_{j=1}^{n}(\bs{x}_{j} - \bar{\bs{x}})(\bs{x}_{j} - \bar{\bs{x}})'
\end{eqnarray}

Then, the generalized likelihood ratio (GLR) is defined by the following equation.
\begin{eqnarray}
\lambda(\bs{x}) & = & \frac{L(\hat{\bs{\mu}}_{1},\hat{\bs{\mu}}_{2},\hat{k},\hat{\bs{\Sigma}})}{L(\hat{\bs{\mu}},\hat{\bs{\Sigma}})}
\end{eqnarray}
Where $-2\log{\lambda(\bs{x})} \sim \chi^{2}(2)$.  The null hypothesis will be rejected if $-2\log{\lambda(\bs{x})} \geq \chi_{1-\alpha}^{2}(2)$.

\section{Constant Variance}
Much like in the previous section, a generalized likelihood test can be used to test for a change in the variance.  In this case we will be testing the following hypothesis.
\begin{equation}
H_{0} : \bs{\Sigma}_{1} = \bs{\Sigma}_{2} = \cdots = \bs{\Sigma}_{k} \neq \bs{\Sigma}_{k+1} = \cdots = \bs{\Sigma}_{n} \qquad H_{a} : \bs{\Sigma}_{1} = \bs{\Sigma}_{2} = \cdots = \bs{\Sigma}_{n} = \bs{\Sigma}
\label{var.increase}
\end{equation}
Likelihood functions for each of the hypotheses will have to be calculated.  To do this, we refer again to equations (\ref{eq.null.llh}) and (\ref{eq.null}).  For the null hypothesis, we have the following likelihood function.
\begin{eqnarray}
L(\bs{\mu},k,\bs{\Sigma}_{1},\bs{\Sigma}_{2}) & = & \prod_{i=1}^{k} f(\bs{x}_{i};\bs{\mu},\bs{\Sigma}_{1}) \prod_{j=k+1}^{n} f(\bs{x}_{j};\bs{\mu},\bs{\Sigma}_{2}) \\
 & = & L(\bs{\mu},\bs{\Sigma}_{1}) \, L(\bs{\mu},\bs{\Sigma}_{2}) 
\end{eqnarray}
Which again reduces to a problem such as exercise (6.11) from our text.  It should be fairly clear that the MLE for $\bs{\mu}$ is $\hat{\bs{\mu}} = \bar{\bs{x}}$, due to the independence of $\bar{\bs{X}}$ and $\bs{S}$.  Then, using the argument in Result 4.10 from the text we can find the MLE's for $\bs{\Sigma}_{1}$ and $\bs{\Sigma}_{2}$, which follow.
\begin{eqnarray}
\hat{\bs{\Sigma}}_{1} & = & \left(\frac{k-1}{k}\right)\sum_{j=1}^{k}(\bs{x}_{j} - \bar{\bs{x}})(\bs{x}_{j} - \bar{\bs{x}})' \\
\hat{\bs{\Sigma}}_{2} & = & \left(\frac{n-k-2}{n-k}\right)\sum_{j=k+1}^{n}(\bs{x}_{j} - \bar{\bs{x}})(\bs{x}_{j} - \bar{\bs{x}})'
\end{eqnarray}
Again, the parameter $k$ cannot be found theoretically, so $\hat{k}$ will be calculated from the data set as the largest value of $L(\hat{\bs{\mu}},k,\hat{\bs{\Sigma}}_{1},\hat{\bs{\Sigma}}_{2})$.

The likelihood function for the alternative hypothesis is simply that found in the previous section defined by equations (\ref{alt}) - (\ref{alt.Sigma}).  The GLR, in this case is defined by the equation
\begin{eqnarray}
\lambda(\bs{x}) & = & \frac{L(\hat{\bs{\mu}},\hat{k},\hat{\bs{\Sigma}}_{1},\hat{\bs{\Sigma}}_{2})}{L(\hat{\bs{\mu}},\hat{\bs{\Sigma}})}
\end{eqnarray}
Where $-2\log{\lambda(\bs{x})} \sim \chi^{2}(2)$ and we reject for large values of $-2\log{\lambda(\bs{x})}$.

\section{Correlation Between Monthly Averages}
Sample correlations between components $X_{i}$ and $X_{k}$ will be calculated using the following formula.
\begin{eqnarray}
r_{ik} & = & \frac{\displaystyle \sum_{j=1}^{n}(x_{ji} - \bar{x}_{i})(x_{jk} - \bar{x}_{k})}{\displaystyle\sqrt{\sum_{j=1}^{n}(x_{ji} - \bar{x}_{i})^{2}}\sqrt{\sum_{j=1}^{n}(x_{jk} - \bar{x}_{k})^{2}}}
\end{eqnarray}
The Pearson correlation test is to test the hypothesis
\begin{equation}
H_{0} : \rho = 0 \textrm{\hspace*{2cm}} H_{a} : \rho \neq 0
\end{equation}
or that the correlation is significant.  The Pearson test uses the following test statistic.
\begin{eqnarray}
t & = & r \sqrt{\frac{n-2}{1-r^{2}}}
\end{eqnarray}
Where $t \sim t_{n-1}$.  This test can be performed simply using \texttt{proc corr} in SAS.  The null hypothesis will be rejected for large $p$-values.

\section{The Normality Assumption}
Johnson and Wichern suggest several ways to assess whether or not the data is actually normally distributed.  The need to verify the normality of the data is not certain because the data is actually the mean vector of the daily temperatures and a central limit argument can be made that the monthly averages should be normally distributed.

Nevertheless, three methods will be used to verify the normality of the data set.  One of the most effective ways is to look at Q-Q plots of each of the components.  If the Q-Q plots turn out to be very close to linear, then normality can be assumed.  Q-Q plots should not show any pattern when varying around the straight line and should not contain an excessive number of points at the tails that stray too far from the fitted line.  Johnson and Wichern also suggest performing a test on the correlation between the data points and the quantile values, this amounts to performing a Shapiro-Wilk test on each of the components.  These tests will be performed.

The text also suggests examining Q-Q plots of the first and last principal components, this seems like a good idea since most of the variability in the data will be contained in the first principal component and almost none in the last.  The implication being that if the first and last principal components are normal, the rest should be well.

\section{Grouping the Months}
Principle component analysis can be used to group the months together.  According to our text, we can determine the principle components by using the eigenvalues and eigenvectors of the sample covariance matrix.  If we have the eigenvalues $\lambda_{1}, \lambda_{2}, \ldots, \lambda_{p}$ with corresponding eigenvectors $\bs{e}_{1}, \bs{e}_{2}, \ldots, \bs{e}_{p}$, then for a random sample $\bs{X} = (X_{1},X_{2},\ldots,X_{p})$, transform $\bs{X}$ to get the new components that follow.
\begin{eqnarray*}
Y_{1} & = & \bs{e}_{1}'\bs{X} \\
Y_{2} & = & \bs{e}_{2}'\bs{X} \\
 & \vdots & \\
Y_{p} & = & \bs{e}_{p}'\bs{X}
\end{eqnarray*}
Call the $Y_{i}$ the principle components and calculate how much of the overall variance is due to each component with the following equation.
\begin{eqnarray}
\left(\begin{array}{c}
\textrm{Proportion of total}\\
\textrm{population variance}\\
\textrm{due to }k \textrm{th principal}\\
\textrm{component}
\end{array}\right)
& = & 
\frac{\lambda_{k}}{\lambda_{1} + \lambda_{2} + \cdots \lambda_{p}} \label{eq.principal}
\end{eqnarray}
Then take the principal components that describe 85\% - 90\% of the variability and use them for all analysis because they are the most important.  To determine this simply add (\ref{eq.principal}) until the sum is larger than 0.85 or 0.90.  It would not be necessary to normalize the data since each column is measured on the same scale.

This analysis can be done easily with \texttt{R} which has functions to calculate eigenvalues and eigenvectors.  There are also several functions available to help in the analysis.

\newpage
\input{bib.tex}
\end{document}
