% \VignetteIndexEntry{Asset selection with Local Search}
% \VignetteKeyword{Local Search}
% \VignetteKeyword{asset allocation}
% \VignetteKeyword{heuristics}
% \VignetteKeyword{optimize}
\documentclass[a4paper]{article}
\usepackage[left=2.5cm,top=2cm, bottom=3cm, right=3.5cm]{geometry}
\usepackage[noae]{Sweave}
\usepackage{mathptmx}
\usepackage{amsmath,amstext}
\usepackage{hyperref}
\usepackage{natbib}
\usepackage{units}
\usepackage{color}
\definecolor{grau2}{rgb}{.2,.2,.2}
\definecolor{grau7}{rgb}{.7,.7,.7}
% define *Sweave* layout
\DefineVerbatimEnvironment{Sinput}{Verbatim}{}
\DefineVerbatimEnvironment{Soutput}{Verbatim}{frame=single,xleftmargin=0em,%
  formatcom=\color{grau2},rulecolor=\color{grau7}}
\DefineVerbatimEnvironment{Scode}{Verbatim}{xleftmargin=2em}
\fvset{listparameters={\setlength{\topsep}{0pt}}}
\renewenvironment{Schunk}{\vspace{\topsep}}{\vspace{\topsep}}
<<echo=false>>=
options(continue = "  ", digits = 5)
@

\begin{document}
{\raggedright{\LARGE Asset selection with Local Search}}\medskip

\noindent Enrico Schumann\\
\noindent \texttt{es@enricoschumann.net}\\
\bigskip


\section{Introduction}

\noindent We provide a code example for a simple
asset-selection model; please see \citet[ch. 12
and~13]{Gilli2011b} for more details.  The code example
here differs slightly from the book's presentation.  To
see the latter, check the script \texttt{exampleLS.R}
(after attaching the package):
<<eval=false>>=
showExample("exampleLS", chapter = "Portfolio")
@
We start by attaching the package and fixing a seed.
<<>>=
require("NMOF")
set.seed(112233)
@

\section{The model}

\noindent We wish to select between $K_{\min}$
and~$K_{\max}$ out of $n_{\mathrm{A}}$ assets such that
an equally-weighted portfolio of these assets has the
lowest-possible variance.  The formal model is:
\begin{align}\label{Eeq:portfolio:modelSelect}
\min_{w}\          & w' \Sigma w &
\intertext{subject to the constraints}
\nonumber w_j      & = \nicefrac{1}{K} \ \text{\quad  for }  j \in J\, , \\[-0.25ex]
\nonumber K_{\min} & \leq K  \leq K_{\max} \, .
\end{align}
The weights are stored in the vector $w$; the symbol
$J$ stands for the set of assets in the portfolio; and
$K = \#\{J\}$ is the cardinality of this set, ie, the
number of assets in the portfolio.


\section{Setting up the algorithm}

We simulate 500~assets: each gets a random volatility
between 20\% and~40\%, and all pairwise correlations
are set to~0.6.
<<>>=
na <- 500L                          ## number of assets
C <- array(0.6, dim = c(na, na))    ## correlation matrix
diag(C) <- 1
minVol <- 0.20; maxVol <- 0.40      ## covariance matrix
Vols <- (maxVol - minVol) * runif(na) + minVol
Sigma <- outer(Vols, Vols) * C
@
The objective function.
<<>>=
OF <- function(x, Data) {
    w <- x/sum(x)
    res <-  crossprod(w[x], Data$Sigma[x, x])
    tcrossprod(w[x], res)
}
@
\ldots or even simpler:
<<>>=
OF2 <- function(x, Data) {
    w <- 1/sum(x)
    sum(w * w * Data$Sigma[x, x])
}
@
The neighbourhood function.
<<>>=
neighbour <- function(xc, Data) {
    xn <- xc
    p <- sample.int(Data$na, Data$nc, replace = FALSE)
    xn[p] <- !xn[p]

    ## reject infeasible solution
    if (sum(xn) > Data$Kmax || sum(xn) < Data$Kmin)
        xc
    else
        xn
}
@
We collect all necessary information in the list
\texttt{Data}: the variance--corvariance matrix
\texttt{Sigma}, the cardinality limits \texttt{Kmin}
and~\texttt{Kmax}, the total number of assets
\texttt{na} (ie, the cardinality of the asset
universe), and the parameter~\texttt{nc}.  This
parameter controls the neighbourhood: it gives the
number of assets that are to be changed when a new
solution is computed.
<<>>=
Data <- list(Sigma = Sigma,
              Kmin = 30L,
              Kmax = 60L,
              na = na,
              nc = 1L)
@

\section{Solving the model}

As an initial solution we use a random portfolio.
<<>>=
card0 <- sample(Data$Kmin:Data$Kmax, 1L, replace = FALSE)
assets <- sample.int(na, card0, replace = FALSE)
x0 <- logical(na)
x0[assets] <- TRUE
@
With this implementation we assume that \texttt{Data\$Kmax >
  Data\$Kmin}. (If \texttt{Data\$Kmax} equals \texttt{Data\$Kmin}, then
\texttt{sample} returns a draw from \texttt{1:Data\$Kmin}.)

We collect all settings for the algorithm in a list \texttt{algo}.
<<>>=
algo <- list(x0 = x0,
              neighbour = neighbour,
              nS = 5000L,
              printDetail = FALSE,
              printBar = FALSE)
@
It remains to run the algorithm.
<<fig = TRUE, height = 3.5>>=
system.time(sol1 <- LSopt(OF, algo, Data))
sqrt(sol1$OFvalue)
par(ylog = TRUE, bty = "n", las = 1, tck = 0.01, mar = c(4,4,1,1))
plot(sqrt(sol1$Fmat[ ,2L]), main = "",
     type = "l", ylab = "portfolio volatility", xlab = "iterations")
@

(Recall that the simulated data had volatilities between 20 and 40\%.)

We can also run the search repeatedly with the same starting value.
<<fig = TRUE, height = 3.5>>=
nRuns <- 3L
allRes <- restartOpt(LSopt, n = nRuns, OF, algo = algo, Data = Data)
allResOF <- numeric(nRuns)
for (i in seq_len(nRuns))
    allResOF[i] <- sqrt(allRes[[i]]$OFvalue)
par(bty = "n", las = 1, tck = 0.01, mar = c(4,4,1,1))
plot(ecdf(allResOF), xlab = "x: Portfolio volatility", pch = 21,
     main = "")
@

(We run \texttt{LSopt} only \Sexpr{nRuns} times to keep the build time
for the vignette acceptable. To get more meaningful results you should
increase \texttt{nRuns}.)

\bibliographystyle{plainnat}
\bibliography{NMOF}
\end{document}