%\VignetteIndexEntry{sdpt3r}
% JSS Article Template
\documentclass{article}


%%% Additional Packages %%%

% Algorithms
% http://tex.stackexchange.com/questions/229355/algorithm-algorithmic-algorithmicx-algorithm2e-algpseudocode-confused
\usepackage{algorithm}% http://ctan.org/pkg/algorithms
\usepackage{algpseudocode}% http://ctan.org/pkg/algorithmicx


% Font stuff
\usepackage[T1]{fontenc}% for correct hyphenation and T1 encoding
%\usepackage{lmodern}% latin modern font
\usepackage[american]{babel}% for American English

% Math fonts/symbols
\usepackage{amsmath}% sophisticated mathematical formulas with amstex (includes \text{})
\usepackage{mathtools}% fix amsmath deficiencies
\usepackage{amssymb}% sophisticated mathematical symbols with amstex
\usepackage{amstext}% ams mathematical fonts
\usepackage{amsfonts}% ams mathematical fonts

% Figures
\usepackage{graphicx}% for including figures

%%% Commands %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% vetors
\newcommand{\ve}[1]{\mathbf{#1}}           % for vetors
\newcommand{\sv}[1]{\boldsymbol{#1}}   % for greek letters
\newcommand{\m}[1]{\mathbf{#1}}               % for matrices
\newcommand{\sm}[1]{\boldsymbol{#1}}   % for greek letters
\newcommand{\tr}[1]{{#1}^{\mkern-1.5mu\mathsf{T}}}              % for transpose
\newcommand{\norm}[1]{||{#1}||}              % for transpose
\newcommand*{\mve}{\operatorname{ve}}
\newcommand*{\trace}{\operatorname{trace}}
\newcommand*{\rank}{\operatorname{rank}}
\newcommand*{\diag}{\operatorname{diag}}
\newcommand*{\vspan}{\operatorname{span}}
\newcommand*{\rowsp}{\operatorname{rowsp}}
\newcommand*{\colsp}{\operatorname{colsp}}
\newcommand*{\svd}{\operatorname{svd}}
\newcommand*{\edm}{\operatorname{edm}}  % euclidean distance matrix (D * D)


% statistical
\newcommand{\widebar}[1]{\overline{#1}}  

% 
% operators
\newcommand{\Had}{\circ}
\DeclareMathOperator*{\lmin}{Minimize}
\DeclareMathOperator*{\argmin}{arg\,min}
\DeclareMathOperator*{\argmax}{arg\,max}
\DeclareMathOperator*{\arginf}{arg\,inf}
\DeclareMathOperator*{\argsup}{arg\,sup}
%\newcommand*{\arginf}{\operatorname*{arginf}}
%\newcommand*{\argsup}{\operatorname*{argsup}}

% Sets
\newcommand*{\intersect}{\cap}
\newcommand*{\union}{\cup}
\let\oldemptyset\emptyset
\let\emptyset\varnothing

% Fields, Reals, etc. etc
\newcommand{\field}[1]{\mathbb{#1}}
\newcommand{\Reals}{\field{R}}
\newcommand{\Integers}{\field{Z}}
\newcommand{\Naturals}{\field{N}}
\newcommand{\Complex}{\field{C}}
\newcommand{\Rationals}{\field{Q}}

% Hyphenation
\hyphenation{Ar-chi-me-dean}

% Editorial
\newcommand*{\TODO}[1]{\textcolor{red}{TODO: #1}}
\newcommand*{\NOTE}[1]{\textcolor{blue}{Note: #1}}
\newtheorem{theorem}{Theorem}
\newtheorem{Proof}{Proof}

\usepackage[total={6.5in,8.75in}, top=0.9in, left=0.7in, right=0.7in, includefoot]{geometry}


% Misc
%\makeatletter
%\newcommand\myisodate{\number\year-\ifcase\month\or 01\or 02\or 03\or 04\or 05\or 06\or 07\or 08\or 09\or 10\or 11\or 12\fi-\ifcase\day\or 01\or 02\or 03\or 04\or 05\or 06\or 07\or 08\or 09\or 10\or 11\or 12\or 13\or 14\or 15\or 16\or 17\or 18\or 19\or 20\or 21\or 22\or 23\or 24\or 25\or 26\or 27\or 28\or 29\or 30\or 31\fi}% create iso date
\makeatother
%\newcommand*{\abstractnoindent}{}% define abstract such that it has no indent
%\let\abstractnoindent\abstract
%\renewcommand*{\abstract}{\let\quotation\quote\let\endquotation\endquote
%  \abstractnoindent}
%\deffootnote[1em]{1em}{1em}{\textsuperscript{\thefootnotemark}}% setting for footnote

\author{Adam Rahman}
\title{The Nearest Correlation Matrix}


%% need no \usepackage{Sweave.sty}

\begin{document}
\SweaveOpts{concordance=FALSE}

\maketitle

First addressed by \cite{higham2002computing} in dealing with correlations between stock prices, difficulty arises when data is not available for all stocks on each day, which is unfortunately a common occurrence. To help address this situation, correlations are calculated for pairs of stocks only when data is available for both stocks on any given day. The resulting correlation matrix is only approximate in that it is not necessarily positive semidefinite.

This problem was cast by \cite{higham2002computing} as

\[
\begin{array}{ll}
\underset{\m{X}}{\text{minimize}} & ||\m{R} - \m{X}||_{F}\\
\text{subject to} & \\
& \begin{array}{rl}
diag(\m{X}) &=~~\ve{1} \\
\m{X} &\in~~\mathcal{S}^{n}
\end{array}
\end{array}
\]

\noindent where $\m{R}$ is the approximate correlation matrix and $||\cdot||_{F}$ denotes the Frobenius norm. Unfortunately, the Frobenius norm in the objective function prevents the problem being formatted as a conic linear optimization problem.

Since the matrix $\m{X}$ is constrained to have unit diagonal and be symmetric, and the matrix $\m{R}$ is an approximate correlation matrix, meaning it will also have unit diagonal and be symmetric, we can re-write the objective function as

\[
||\m{R}-\m{X}||_{F} = 2*||svec(\m{R}) - svec(\m{X})|| = 2*||\ve{e}||
\]

Now, introduce a variable $e_{0}$ such that $e_{0} \geq ||\ve{e}||$, and define $\ve{e}^{\ast} = [e_{0};~\ve{e}]$. The vector $\ve{e}^{\ast}$ is now restricted to be in the quadratic cone $\mathcal{Q}^{n(n+1)/2+1}$. This work leads to the formulation of \cite{toh1999sdpt3}

\[
\begin{array}{ll}
\underset{\ve{e}^{\ast},~\m{X}}{\text{minimize}} & e_{0} \\
\text{subject to} & \\
& 
\begin{array}{rl}
svec(\m{R}) - svec(\m{X}) &=~~ [\ve{0}, \m{I}_{n(n+1)/2}]~\ve{e}^{\ast} \\
diag(\m{X}) &=~~\ve{1} \\
\m{X} &\in~~\mathcal{S}^{n} \\
\ve{e}^{\ast} &\in~~\mathcal{Q}^{n(n+1)/2+1}
\end{array}
\end{array}
\]

Here, $[\m{X},\m{Y}]$ denotes column binding of the two matrices $\m{X}_{n\times p}$ and $\m{Y}_{n \times m}$ to form a matrix of size $n \times (p+m)$. By minimizing $e_{0}$, we indirectly minimize $\ve{e} = svec(\m{R}) - svec(\m{X})$, since recall we have $e_{0} \geq ||\ve{e}||$, which is the goal of the original objective function.

To see this as a conic linear optimization problem, notice that $e_{0}$ can be written as $\langle \m{C}^{q},~\m{X}^{q} \rangle$ by letting $\m{C}^{q} = [1;\ve{0}_{n(n+1)/2}]$ and $\m{X}^{q} = \ve{e}^{\ast}$. Since the matrix $\m{X}$ (i.e. $\m{X}^{s}$) does not appear in the objective function, the matrix $\m{C}^{s}$ is an $n \times n$ matrix of zeros.

Re-writing the first constraint as

\[
svec(\m{X}) + [\ve{0}, \m{I}_{n(n+1)/2}]~\ve{e}^{\ast} =~~  svec(\m{R})
\]

\noindent we can easily define the constraint matrices and right hand side of the first constraint as

\[
\begin{array}{rl}
\m{A}_{1}^{s} &= \m{I}_{n(n+1)/2} \\
\m{A}_{1}^{q} &= [\ve{0}, \m{I}_{n(n+1)/2}] \\
\ve{b}_{1} &= svec(\m{R})
\end{array}
\]

The second constraint is identical to the constraint from the Max-Cut problem, where each diagonal element of $\m{X}$ is constrained to be equal to 1. Define $\ve{b}_{2} = \ve{1}$, and for the $k^{th}$ diagonal element of $\m{X}$, define the matrix $\m{A}_{k}$ as

\[
\m{A}_{k} = [a_{ij}] = 
\begin{cases}
1, & i = j = k \\
0, & \text{otherwise}
\end{cases}
\]

\noindent yielding $\langle \m{A}_{k},\m{X} \rangle = x_{kk}$. To write this as $\tr{(\m{A}_{2}^{s})}\m{X}^{s}$, define 

\[
\m{A}_{2}^{s} = [svec(\m{A}_{1}),...,svec(\m{A}_{n})]
\]

Since $\ve{e}^{\ast}$ does not appear in the second constraint, $\m{A}_{2}^{q} = \ve{0}_{n(n+1)/2+1}$. 

The final step is to combine the individual constraint matrices from each constraint to form one constraint matrix for each variable, which is done by defining $\m{A}^{s} = [\m{A}_{1}^{s},~\m{A}_{2}^{s}]$, $\m{A}^{q} = [\m{A}_{1}^{q},~\m{A}_{2}^{q}]$. We also concatenate both right hand side vectors to form a single vector by defining $\ve{b} = [\ve{b}_{1};~\ve{b}_{2}]$. Here, the notation $[\m{X};\m{Y}]$ is used to denote two matrices $\m{X}_{p \times m}$ and $\m{Y}_{q \times m}$ bound vertically to form a matrix of size $(p+q) \times m$. With this, the nearest correlation matrix problem is written as a conic linear optimization.

To solve this problem using \verb!sqlp!, we first define \verb!blk!. There are two optimization variables in the formulation of the nearest correlation matrix - $\m{X}$ is an $n \times n$ matrix constrained to be in a semidefinite cone, and $\ve{y}$ is an $n(n+1)/2 + 1$ length vector constrained to be in a quadratic cone, so 

\begin{verbatim}
R> blk <- matrix(list(),nrow=2,ncol=2)

R> blk[[1,1]] <- "s"
R> blk[[1,2]] <- n

R> blk[[2,1]] <- "q"
R> blk[[2,2]] <- n * (n + 1) / 2 + 1
\end{verbatim}

Note that $\m{X}$ does not appear in the objective function, so the \verb!C! entry corresponding to the block variable $\m{X}$ is an $n \times n$ matrix of zeros, which defines \verb!C! as

\begin{verbatim}
R> C <- matrix(list(),nrow=2,ncol=1)
R> C[[1,1]] <- matrix(0,nrow=n,ncol=n)
R> C[[2,1]] <- rbind(1, matrix(0,nrow=n2,ncol=1))
\end{verbatim}

Next comes the constraint matrix for $\m{X}$

\begin{verbatim}
R> At <- matrix(list(),nrow=2,ncol=1)
R>
R> #Constraint Matrix for Upper Triangular Elements of X
R> A1s <- diag(1,nrow=n2,ncol=n2)
R>
R> #Construct Ak matrices
R> Aks <- matrix(list(),nrow=1,ncol=n)
R> for(k in 1:n){
R>   Aks[[k]] <- matrix(0,nrow=n,ncol=n)
R>   diag(Aks[[k]])[k] <- 1
R> }
R>
R> A2s <- svec(blk[1,],Aks)
R> 
R> #Combined Constraint Matrix for X
R> At[[1,1]] <- cbind(A1s,A2s)
\end{verbatim}

\noindent then the constraint matrix for $\ve{e}^{\ast}$.

\begin{verbatim}
R> A1q<- matrix(0,nrow=n,ncol=n2+1)
R>
R> A2q1 <- matrix(0,nrow=n2,ncol=1)
R> A2q2 <- diag(1,nrow=n2,ncol=n2)
R> A2q <- cbind(A211, A212)
R>
R> At[[2,1]] <- rbind(A1q, A2q)
\end{verbatim}

\noindent and the right hand side vector \verb!b!

\begin{verbatim}
R> b <- rbind(svec(blk[1,],R),matrix(1,n,1))
\end{verbatim}

The nearest correlation matrix problem is now solved by

\begin{verbatim}
R> sqlp(blk, At, C, b)
\end{verbatim}

To demonstrate the nearest correlation matrix problem, we will modify an existing correlation matrix by exploring the effect of changing the sign of just one of the pairwise correlations. In the context of stock correlations, we make use of tools available in the \verb!R! package \verb!quantmod! (\cite{quantmod}) to access stock data from five tech firms (Microsoft, Apple, Amazon, Alphabet/Google, and IBM) beginning in 2007.

\begin{verbatim}
R> library(quantmod)

R> getSymbols(c("MSFT", "AAPL", "AMZN", "GOOGL", "IBM"))
R> stock.close <- as.xts(merge(MSFT, AAPL, AMZN, GOOGL,IBM))[,c(4,10,16,22,28)]
\end{verbatim}

The correlation matrix for these five stocks is

\begin{verbatim}
R> stock.corr <- cor(stock.close)
R> stock.corr


            MSFT.Close AAPL.Close AMZN.Close GOOGL.Close IBM.Close
MSFT.Close   1.0000000 -0.2990463  0.9301085   0.5480033 0.2825698
AAPL.Close  -0.2990463  1.0000000 -0.1514348   0.3908624 0.6887127
AMZN.Close   0.9301085 -0.1514348  1.0000000   0.6228299 0.3870390
GOOGL.Close  0.5480033  0.3908624  0.6228299   1.0000000 0.5885146
IBM.Close    0.2825698  0.6887127  0.3870390   0.5885146 1.0000000
\end{verbatim}

Next, consider the effect of having a positive correlation between Microsoft and Apple

\begin{verbatim}
R> stock.corr[1,2] <- -1 * stock.corr[1,2]
R> stock.corr[2,1] <- stock.corr[1,2]
R> stock.corr

            MSFT.Close AAPL.Close AMZN.Close GOOGL.Close IBM.Close
MSFT.Close   1.0000000  0.2990463  0.9301085   0.5480033 0.2825698
AAPL.Close   0.2990463  1.0000000 -0.1514348   0.3908624 0.6887127
AMZN.Close   0.9301085 -0.1514348  1.0000000   0.6228299 0.3870390
GOOGL.Close  0.5480033  0.3908624  0.6228299   1.0000000 0.5885146
IBM.Close    0.2825698  0.6887127  0.3870390   0.5885146 1.0000000
\end{verbatim}

Unfortunately, this correlation matrix is not positive semidefinite

\begin{verbatim}
R> eigen(stock.corr)$values

[1]  2.8850790  1.4306393  0.4902211  0.3294150 -0.1353544
\end{verbatim}

Given the approximate correlation matrix \verb!stock.corr!, the built-in function \verb!nearcorr! provides the input variables required for the problem to be solved using \verb!sqlp!

\begin{verbatim}
R> out <- nearcorr(stock.corr)
R> blk <- out$blk
R> At <- out$At
R> C <- out$C
R> b <- out$b

R> foo <- sqlp(blk,At,C,b)
\end{verbatim}

Since this is a minimization problem, and thus a primal formulation of the SQLP, the output \verb!X! from \verb!sqlp! will provide the optimal solution to the problem - that is, \verb!X! will be the nearest correlation matrix to \verb!stock.corr!.

\begin{verbatim}
foo$X

          [,1]        [,2]        [,3]      [,4]      [,5]
[1,] 1.0000000  0.25388359  0.86150833 0.5600734 0.3126420
[2,] 0.2538836  1.00000000 -0.09611382 0.3808981 0.6643566
[3,] 0.8615083 -0.09611382  1.00000000 0.6115212 0.3480430
[4,] 0.5600734  0.38089811  0.61152116 1.0000000 0.5935021
[5,] 0.3126420  0.66435657  0.34804303 0.5935021 1.0000000
\end{verbatim}

The matrix above is symmetric with unit diagonal and all entries in $[-1,1]$. By checking the eigenvalues,

\begin{verbatim}
eigen(X)

$values

[1] 2.846016e+00 1.384062e+00 4.570408e-01 3.128807e-01 9.680507e-11
\end{verbatim}

\noindent we can see that \verb!X! is indeed a correlation matrix.

\bibliographystyle{plain}
\bibliography{sdpt3r}
\end{document}