The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
34,139 characters
Implementing Convex Optimization in R: Two Econometric Examples
\title{Implementing Convex Optimization in R:\\ Two Econometric Examples}
\author{Zhan Gao \\ \small University of Southern California \vspace{5pt}\\ Zhentao Shi\\ \small The Chinese University of Hong Kong}
\date{}
\maketitle
\vspace{-.5cm}
\begin{abstract}
Economists specify high-dimensional models to address heterogeneity
in empirical studies with complex big data. Estimation of these models
calls for optimization techniques to handle a large number of parameters.
Convex problems can be effectively executed in modern statistical
programming languages. We complement Koenker and Mizera (2014)'s work
on numerical implementation of convex optimization, with focus on
high-dimensional econometric estimators. Combining \texttt{R} and
the convex solver \texttt{MOSEK} achieves faster speed and equivalent
accuracy, demonstrated by examples from Su, Shi, and Phillips (2016)
and Shi (2016). Robust performance of convex optimization is witnessed
cross platforms. The convenience and reliability of convex optimization
in \texttt{R} make it easy to turn new ideas into prototypes.
\end{abstract}
\vspace{1cm}
\noindent Key words: big data, convex optimization, high-dimensional
model, numerical solver
\noindent JEL code: C13, C55, C61, C87
\vspace{1cm}
\small \noindent Zhan Gao: \texttt{[email removed]}. Address: Department
of Economics, University of Southern California, 3620 South Vermont
Ave. Kaprielian Hall, 300 Los Angeles, CA 90089-0253, USA. Zhentao
Shi (corresponding author): \texttt{[email removed]}. Address:
Department of Economics, 912 Esther Lee Building, the Chinese University
of Hong Kong, Sha Tin, New Territories, Hong Kong SAR, China. Tel:
(852) 3943-1432. Fax (852) 2603-5805. We thank Roger Koenker for inspiration
and hospitality during the second author's visit to University of
Illinois.
\newpage{}
\normalsize
\section{Introduction}
Equipped with tremendous growth of computing power over the last few
decades, econometricians endeavor to tackle high-dimensional real
world problems that we could hardly have imagined before. Along with
the development of modern asymptotic theory, computation has gradually
ascended onto the central stage. Today, discussion of numerical algorithms
is essential for new econometric procedures.
Optimization is at the heart of estimation, and convex optimization
is the best understood category. Convex problems are ubiquitous in
econometric textbooks. The least square problem is convex, and the
classical normal regression is also convex after straightforward reparametrization.
Given a linear single-index form, the Logit or Probit binary regression,
the Poisson regression and the regressions with a censored or truncated
normal distributions are all convex. Another prominent example is
the quantile regression \citep{koenker1978regression}, motivated
from its robustness to non-Gaussian errors and outlier contamination.
With the advent of big data, practitioners attempt to build general
models that involve hundreds or even more parameters in the hope to
capture complex heterogeneity in empirical economic studies. Convex
optimization techniques lay out the foundation of estimating these
high-dimensional models. Recent years witnesses \citet{bajari2015machine},
\citet{gu2017empirical} and \citet{doudchenko2016balancing}, to
name a few, exploring new territories by taking advantage of convexity.
To facilitate practical implementation, \citet{cvxinr} summarize
the operation in \texttt{R} by \texttt{MOSEK} via \texttt{Rmosek}
to solve linear programming, conic quadratic programming, quadratic
programming, etc. \texttt{R} is open-source software, \texttt{MOSEK}
is a proprietary convex optimization solver but offers free academic
license, and \texttt{Rmosek} is the \texttt{R} interface that communicates
with \texttt{MOSEK}. \texttt{MOSEK} specializes in convex problems
with reliable performance, and is competitive in high-dimensional
problems.
This paper complements \citet{cvxinr}'s work. We replicate by \texttt{Rmosek}
two examples of high-dimensional estimators, namely \citet{su2016identifying}'s
classifier-Lasso (C-Lasso) and \citet{REL}'s relaxed empirical likelihood
(REL). In addition, we replicate by \texttt{Rmosek} an application
of C-Lasso that re-examines China's GDP growth rate \citep{chinagdp}.
These exercises highlight two points. Firstly, the \texttt{R} environment
is robust in numerical accuracy for high-dimensional convex optimization
and \texttt{Rmosek} takes the lead in computational speed. Second,
we showcase the ease of creating new econometric estimators\textemdash often
no more than a few lines of code\textemdash by the code snippets in
the Appendix. Such convenience lowers the cost of turning an idea
into a prototype, and enables researchers to glean valuable insights
about their archetypes by experimenting new possibilities. All code
in this note is hosted at \url{https://github.com/zhentaoshi/convex_prog_in_econometrics}.
\section{Classifier-Lasso \label{sec:Classifier-Lasso}}
It is common practice to assume in linear fixed-effect panel data
models that the cross-sectional units are heterogeneous in terms of
the time-invariant individual intercept, while they all share the
same slope coefficient. This pooling assumption can be tested and
is often rejected in real-world applications. In recent years panel
data group structure has been attracting attention. \citet{bonhomme2015grouped}
allow group structure in the intercept and use the $k$-means algorithm
for classification. When the slope coefficients exhibit group structure,
\citet{su2016identifying} propose Classifier-Lasso (C-Lasso) to identify
the latent group pattern.
We illustrate the penalized least square (PLS), a simple special case
of C-Lasso.\footnote{The profile log-likelihood function $Q_{1,nT}\left(\beta\right)=\sum_{i=1}^{n}\sum_{t=1}^{T}\psi\left(w_{it},\beta_{i},\hat{\mu}_{i}\left(\beta_{i}\right)\right)$
for nonlinear models can be reformulated into a separable form, while
penalized GMM (PGMM) can be handled under the same optimization framework
as PLS. They are discussed in Appendix \ref{subsec:Nonlinear-Lasso}
and Appendix \ref{subsec:Penalized-GMM}, respectively.} Given a tuning parameter $\lambda$ and the number of groups $K$,
PLS is defined as the solution to
\[
\min_{\boldsymbol{\beta},\left(\alpha_{k}\right)_{k=1}^{K}}\,\frac{1}{nT}\sum_{i=1}^{n}\sum_{t=1}^{T}\left(y_{it}-x_{it}^{\prime}\beta_{i}\right)^{2}+\frac{\lambda}{n}\sum_{i=1}^{n}\prod_{k=1}^{K}\left\Vert \beta_{i}-\alpha_{k}\right\Vert _{2}
\]
where $\boldsymbol{\beta}=\left(\beta_{i}\right)_{i=1}^{n}$. The
additive-multiplicative penalty pushes the individual slope coefficients
$\beta_{i}$ in the same group toward a common coefficient $\alpha_{k}$.
This is not a convex problem, but the optimization with the additive-multiplicative
penalty can be approximated by an iterative algorithm, as is explained
in the Supplement of \citet[Section S3.1]{su2016identifying}. Procedures
based on such an iteration have been successfully applied to \citet{su2017interactive},
\citet{su2017numberofgroup} and \citet{su2017sieve}. The iterative
algorithm initiates at the within-group estimator, which is consistent
when $T$ is large. In the $k$-th sub-step of the $r$-th iteration,
$(\boldsymbol{\beta},\alpha_{\tilde{k}})$ is chosen to minimize
\begin{equation}
\min_{\boldsymbol{\beta},\alpha_{\tilde{k}}}\,\frac{1}{nT}\sum_{i=1}^{n}\sum_{t=1}^{T}\left(y_{it}-x_{it}^{\prime}\beta_{i}\right)^{2}+\frac{\lambda}{n}\sum_{i=1}^{n}\left\Vert \beta_{i}-\alpha_{\tilde{k}}\right\Vert _{2}\gamma_{i}\label{eq:iter}
\end{equation}
where $\gamma_{i}=\prod_{k=1}^{\tilde{k}-1}\Vert\hat{\beta}_{i}^{(r,k)}-\hat{\alpha}_{k}^{(r)}\Vert_{2}\cdot\prod_{k=\tilde{k}+1}^{K}\Vert\hat{\beta}_{i}^{(r-1,k)}-\hat{\alpha}_{k}^{(r-1)}\Vert_{2}.$
The iteration proceeds until the $K$-convex problem numerically converges.
Given the multiplier $\gamma_{i}$, the above optimization problem
is convex in $(\boldsymbol{\beta},\alpha_{\tilde{k}})$ and the structure
is very close to Lasso. Though the \texttt{R} packages \texttt{lars}
and\texttt{ glmnet} packages can carry out the standard Lasso, however,
it is not straightforward how to modify these functions to accommodate
(\ref{eq:iter}), where $\alpha_{\tilde{k}}$ is also an unknown parameter
to be optimized. A quick review of \citet{cvxinr} approach to Lasso
will be helpful.
\bigskip
\begin{example}[Lasso]
The standard Lasso problem is
\[
\min_{\beta}\,\frac{1}{n}\left\Vert y-X\beta\right\Vert _{2}^{2}+\lambda\left\Vert \beta\right\Vert _{1}
\]
where $y\in\mathbb{R}^{n}$ and $X\in\mathbb{R}^{n\times p}$ are
observed data, $\lambda$ is the tuning parameter and $\beta\in\mathbb{R}^{p}$
is the parameter of interests. However, \texttt{Rmosek} does not accept
the $l_{1}$ norm. To overcome the difficulty, \citet{cvxinr} introduce
new parameters to transform the $l_{1}$-penalized problem into a
conic optimization that \texttt{Rmosek} recognizes. We first deal
with $\left\Vert \beta\right\Vert _{1}$. The $p\times1$ vector $\beta$
can be decomposed into a positive part $\beta^{+}=\left(\max\left\{ 0,\beta_{j}\right\} \right)_{j=1}^{p}$
and a negative part $\beta^{-}=\left(\max\left\{ 0,-\beta_{j}\right\} \right)_{j=1}^{p}$,
so that $\beta=\beta^{+}-\beta^{-}$ and $\left\Vert \beta\right\Vert _{1}=e^{\prime}\beta^{+}+e^{\prime}\beta^{-}$,
where $e$ is the $p\times1$ vector with all elements equal to $1$.
Next, we transform the $l_{2}$-norm $\left\Vert y-X\beta\right\Vert _{2}^{2}$
to a second-order conic constraint. Consider a minimization problem
with $\left\Vert v\right\Vert _{2}^{2}$ in the objective function.
We can use a new parameter $t$ to replace it and add a conic constraint
$\left\Vert v\right\Vert _{2}^{2}\leq t$, which is equivalent to
$\left\Vert \left(v,\frac{t-1}{2}\right)\right\Vert _{2}\leq\frac{t+1}{2}$.
Thus we obtain a standard conic constraint $\left\Vert \left(v,s\right)\right\Vert _{2}\leq r$,
where $s=\frac{t-1}{2}$ and $r=\frac{t+1}{2}$. We rewrite the Lasso
problem as
\begin{gather*}
\min_{\theta}\,\lambda\left(e'\beta^{+}+e'\beta^{-}\right)+\frac{t}{n}\\
\text{s.t.\ \ }v=y-X\left(\beta^{+}-\beta^{-}\right),\,\left\Vert \left(v,s\right)\right\Vert _{2}\leq r,\,s=\frac{t-1}{2},\,r=\frac{t+1}{2}
\end{gather*}
where $\theta=\left(\beta^{+},\beta^{-},v,t,s,r\right)$. This problem
is of the standard form of second-order conic programming and hence
can be executed in \texttt{Rmosek}.
\bigskip
\end{example}
Applying the techniques in the Lasso formulation, we can transform
the $l_{2}$-norm terms in (\ref{eq:iter}) and formulate the problem
into a conic programming:
\begin{gather*}
\min_{\alpha_{\tilde{k}},\theta}\,\sum_{i=1}^{n}\left(\left(\frac{1}{nT}\right)t_{i}+\left(\frac{\lambda}{n}\gamma_{i}\right)w_{i}\right)\\
\text{s.t.}\;\;x_{i}\beta_{i}+\nu_{i}=y_{i},\,\beta_{i}-\mu_{i}-\alpha_{\tilde{k}}=0,\,s_{i}-\frac{1}{2}t_{i}=-\frac{1}{2},\:r_{i}-\frac{1}{2}t_{i}=\frac{1}{2},\\
\left\Vert \left(\nu_{i},s_{i}\right)\right\Vert _{2}\leq r_{i},\,\left\Vert \mu_{i}\right\Vert _{2}\leq w_{i},\,t_{i}\geq0,\,\text{for all}\,i=1,2,\cdots,n
\end{gather*}
where $\theta=\left\{ \beta_{i},\nu_{i},\mu_{i},s_{i},r_{i},t_{i},w_{i}\right\} _{i=1}^{n}$.
The convexity is manifest when we write the problem in matrix form,
as is displayed in Appendix \ref{sec:Code-Snippets}.
\subsection{Replication}
We replicate\textcolor{black}{{} }the simulation studies in \citet[Section 4]{su2016identifying}
in \texttt{R} via \texttt{Rmosek} and compare the performance of different
numerical optimization approaches. \citet{su2016identifying} conduct
their numerical work in \texttt{MATLAB} via \texttt{CVX} \citep{cvx}.
\texttt{CVX} is a \texttt{MATLAB} add-on package for \emph{disciplined
convex optimization} \citep[DCP]{grant2006disciplined}. It provides
an interface to communicate with commercial or open-source solvers.
In the \texttt{R }\textcolor{black}{environment, the de facto solver
is }\texttt{\textcolor{black}{optimx}}\textcolor{black}{{} \citep{optimx};
another option is the interface }\texttt{\textcolor{black}{nloptr}}\textcolor{black}{{}
\citep{nloptr} that hooks optimization solver }\texttt{\textcolor{black}{NLopt}}\textcolor{black}{{}
\citep{nlopt}. They are general-purpose optimization solvers not
tailored for convexity. Most recently, \citet{CVXR} are actively
developing }\texttt{\textcolor{black}{CVXR}}\textcolor{black}{, }\texttt{\textcolor{black}{CVX}}\textcolor{black}{'s
counterpart in }\texttt{\textcolor{black}{R}}\textcolor{black}{. At
this stage, it is integrated with the open-source solver }\texttt{\textcolor{black}{ECOS}}\textcolor{black}{{}
\citep{ECOS}. We also consider the counterpart of }\texttt{\textcolor{black}{CVX}}\textcolor{black}{{}
in }\texttt{\textcolor{black}{Python}}\textcolor{black}{{} environment,
}\texttt{\textcolor{black}{CVXPY}}\textcolor{black}{{} \citep{cvxpy},
to verify the stability of the algorithm across platforms.}\footnote{\textcolor{black}{In the latest version (Version 0.99), }\texttt{\textcolor{black}{CVXR}}\textcolor{black}{{}
supports }\texttt{\textcolor{black}{MOSEK}}\textcolor{black}{{} by sending
the problem to }\texttt{\textcolor{black}{MOSEK}}\textcolor{black}{{}
in the }\texttt{\textcolor{black}{Python}}\textcolor{black}{{} environment.
In our experiment, large-scale problems like the C-Lasso cause errors
in the communication between }\texttt{\textcolor{black}{R}}\textcolor{black}{{}
and }\texttt{\textcolor{black}{Python}}\textcolor{black}{. In addition,
}\texttt{\textcolor{black}{CVXR}}\textcolor{black}{{} with }\texttt{\textcolor{black}{MOSEK}}\textcolor{black}{{}
currently cannot incorporate problems with nonlinear objective functions
and hence cannot be used for REL in Section \ref{sec:Relaxed-Empirical-Likelihood}. }}
\begin{table}[h]
\caption{\label{tab:classopls}Classification and Point Estimation of $\alpha_{1}$:
Replication of \citet[DGP 1]{su2016identifying}}
\medskip{}
\centering{}
\begin{tabular}{rrrrrrr}
\toprule
$\left(n,T\right)$ &
$\left(100,15\right)$ &
$\left(100,25\right)$ &
$\left(100,50\right)$ &
$\left(200,15\right)$ &
$\left(200,25\right)$ &
$\left(200,50\right)$\tabularnewline
\midrule
\multicolumn{7}{c}{RMSE}\tabularnewline
\texttt{Rmosek} &
0.0762 &
0.0386 &
0.0247 &
0.0428 &
0.0278 &
0.0174\tabularnewline
\texttt{CVXR} &
0.0762 &
0.0386 &
0.0247 &
0.0427 &
0.0278 &
0.0174\tabularnewline
\texttt{CVX} &
0.0767 &
0.0399 &
0.0253 &
0.0443 &
0.0286 &
0.0179\tabularnewline
\texttt{CVXPY} &
0.0741 &
0.0394 &
0.0253 &
0.0424 &
0.0271 &
0.0173\tabularnewline
\midrule
\multicolumn{7}{c}{Correct Ratio}\tabularnewline
\texttt{Rmosek} &
0.8987 &
0.9645 &
0.9965 &
0.9019 &
0.9668 &
0.9969\tabularnewline
\texttt{CVXR} &
0.8986 &
0.9645 &
0.9965 &
0.9020 &
0.9668 &
0.9969\tabularnewline
\texttt{CVX} &
0.8991 &
0.9647 &
0.9965 &
0.9026 &
0.9667 &
0.9968\tabularnewline
\texttt{CVXPY} &
0.8988 &
0.9644 &
0.9965 &
0.9021 &
0.9667 &
0.9969\tabularnewline
\midrule
\multicolumn{7}{c}{Running Time (in minute) }\tabularnewline
\texttt{Rmosek} &
18.08 &
10.58 &
8.42 &
24.94 &
15.87 &
13.57\tabularnewline
\texttt{CVXR} &
73.85 &
40.75 &
32.51 &
77.86 &
47.26 &
37.03\tabularnewline
\texttt{CVX} &
94.57 &
50.25 &
34.98 &
90.14 &
54.29 &
41.27\tabularnewline
\texttt{CVXPY} &
27.91 &
18.17 &
26.32 &
33.64 &
26.70 &
29.09\tabularnewline
\bottomrule
\end{tabular}
\end{table}
We follow DGP 1 in \citet[Section 4]{su2016identifying} as a benchmark.
Table \ref{tab:classopls} reports under various combinations of the
cross sectional units $n$ and the time length $T$, the root-mean-square
error (RMSE) of $\widehat{\alpha}_{1}$ and the probability of correct
group classification (correct ratio). The DGP, simulation settings
and the indicators are relegated to Appendix \ref{subsec:C-Lasso DGP}
to save space.
Within the \texttt{R} environment, the numerical results of estimation
error and classification correct ratio by \texttt{Rmosek} are almost
identical to \texttt{CVXR} up to rounding errors.\footnote{The de facto solver \texttt{optimx} breaks down when solving such
high dimensional problems. \texttt{nloptr} takes more than a few hours
to finish one estimation, which makes the full-scale simulation exercise
computational infeasible. In addition, \texttt{nloptr} fails to attain
an accurate solution in most cases according to our experiments.} We also implement the simulation in \texttt{MATLAB} via \texttt{CVX}
and in \texttt{Python} via \texttt{CVXPY}, the results are largely
similar, which demonstrates the robustness of the numerical performance
of C-Lasso across different computing platforms.
Practitioners may need to try out different specifications for robustness
check in real applications. Without fast optimization solvers, computational
cost can become a bottleneck. On the same computing platform of Intel(R)
Core(TM) i7-8750H CPU @ 2.20GHz, each case is executed in a single
thread and we record the running time in the lower panel in Table
\ref{tab:classopls}. \texttt{Rmosek} significantly outperforms all
alternatives. \texttt{CVX} in \texttt{MATLAB} is about $3$ to $5.2$
times slower and \texttt{CVXPY} is about $1.4$ to $3.1$ times slower
than \texttt{Rmosek}. Although \texttt{CVX} and \texttt{CVXPY} are
also powered by \texttt{MOSEK}, the DCP system takes time to check
the convexity of the input problem and automate the formulation. For
similar reasons, \texttt{CVXR} is $2.7$ to $4$ times slower than
\texttt{Rmosek}. According to \citet[Section 4.1]{CVXR}, we can skip
the DCP formulation steps with \texttt{CVXR} and the advantage of
\texttt{Rmosek} becomes around $1.5$ to $1.8$, which illustrates
the advantages of \texttt{MOSEK} over the open-source solver \texttt{ECOS}.
In summary, DCP is useful when we are uncertain about the convexity
and solvability of a problem. However, for problems that are mathematically
verified to be convex, directly calling \texttt{MOSEK} saves much
computational time.
\subsection{Empirical Application\label{subsec:Empirical-Application}}
While China is now the second largest economy in the world in terms
of aggregate GDP, the accuracy of its reported national income accounting
has been a topic of constant debate over the years. Most recently,
\citet{chinagdp} utilize local economic indicators that are directly
associated with economic activities to estimate China's local and
aggregate GDP in order to assess the quality of these numbers. Different
regions of this continent-size country are growing at varying pace,
thereby resulting in tremendous heterogeneity among its provinces.
To control the hidden heterogeneity, \citet{chinagdp} specify a linear
fixed effect model with latent group structure
\[
y_{it}=x_{it}^{\prime}\beta_{i}+v_{i}+\varepsilon_{it},
\]
where $y_{it}$ is the logarithm of GDP for province $i$ at year
$t$, $x_{it}$ includes local economic indicators of interests, $v_{i}$
characterizes the fixed effect of province $i$, and $\varepsilon_{it}$
is the idiosyncratic error. The heterogeneous slope coefficients $\beta_{i}$
captures latent group structures across regions to be determined by
C-Lasso. The data span from year 2000 to 2007, i.e. $T=8$. Five indicators
are employed as regressors to control observable heterogeneity, namely
\emph{satellite night lights}, \emph{national tax revenue}, \emph{exports},
\emph{imports}, and \emph{electricity consumption}. These indicators
are less susceptible to local officials' manipulation and thus more
robustly reflect economic activities for real businesses. Two alternative
specifications, one with no \emph{satellite night lights} and the
other with neither \emph{satellite night lights} nor \emph{national
tax revenue}, are also considered.
In our implementation, the number of groups $K$ and tuning parameter
$\lambda$ are determined by the information criterion proposed in
\citet[Section 2.5]{su2016identifying}. Trials with different specifications
and tuning parameters can be time-consuming, particularly for our
panel data of a short $T$ as the algorithm has to iterate many times
until numerical convergence. \citet{chinagdp} estimate the model
in \texttt{MATLAB} via \texttt{CVX} and the classification results
are reported in \citet[Table A11]{chinagdp}. We replicate the classification
results by \texttt{Rmosek} and compare the accuracy and speed to \texttt{CVX}.
\begin{table}
\caption{\label{tab:chinagdp_result_table}Classification Results: Replication
of \citet[Table A11]{chinagdp}}
\medskip{}
\begin{centering}
\begin{tabular}{llllllll}
\toprule
\multicolumn{2}{c}{All 5 Indicators} &
&
\multicolumn{2}{c}{Without Light} &
&
\multicolumn{2}{c}{Without Light and Tax}\tabularnewline
Group 1 &
Group 2 &
&
Group 1 &
Group 2 &
&
Group 1 &
Group 2\tabularnewline
\midrule
{\small{}Beijing} &
{\small{}Tianjin} &
&
{\small{}Beijing} &
{\small{}Tianjin} &
&
{\small{}Beijing} &
{\small{}Tianjin}\tabularnewline
{\small{}Hebei} &
{\small{}Jilin} &
&
{\small{}Inner Mongolia} &
{\small{}Hebei} &
&
{\small{}Hebei} &
{\small{}Liaoning}\tabularnewline
{\small{}Shanxi} &
{\small{}Heilongjiang} &
&
{\small{}Liaoning} &
{\small{}Shanxi} &
&
{\small{}Shanxi} &
{\small{}Shanghai}\tabularnewline
{\small{}Inner Mongolia} &
\textbf{\small{}Liaoning} &
&
{\small{}Jilin} &
{\small{}Heilongjiang} &
&
{\small{}Inner Mongolia} &
{\small{}Zhejiang}\tabularnewline
{\small{}Shanhai} &
{\small{}Jiangxi} &
&
{\small{}Shanghai} &
{\small{}Zhejiang} &
&
{\small{}Jilin} &
{\small{}Shandong}\tabularnewline
{\small{}Jiangsu} &
{\small{}Henan} &
&
{\small{}Jiangsu} &
{\small{}Jiangxi} &
&
{\small{}Heilongjiang} &
{\small{}Henan}\tabularnewline
{\small{}Anhui} &
{\small{}Hunan} &
&
{\small{}Anhui} &
{\small{}Shandong} &
&
{\small{}Jiangsu} &
{\small{}Hunan}\tabularnewline
{\small{}Fujing} &
{\small{}Guangdong} &
&
{\small{}Fujian} &
{\small{}Henan} &
&
{\small{}Anhui} &
{\small{}Guangdong}\tabularnewline
{\small{}Hubei} &
\textbf{\small{}Zhejiang} &
&
{\small{}Hubei} &
{\small{}Hunan} &
&
{\small{}Fujian} &
{\small{}Chongqing}\tabularnewline
{\small{}Hainan} &
{\small{}Guangxi} &
&
{\small{}Hainan} &
{\small{}Guangdong} &
&
{\small{}Jiangxi} &
{\small{}Guizhou}\tabularnewline
{\small{}Qinghai} &
{\small{}Chongqing} &
&
{\small{}Qinghai} &
{\small{}Guangxi} &
&
{\small{}Hubei} &
{\small{}Yunnan}\tabularnewline
{\small{}Xinjiang} &
{\small{}Sichuan} &
&
&
{\small{}Chongqing} &
&
{\small{}Guangxi} &
{\small{}Shaanxi}\tabularnewline
&
\textbf{\small{}Shandong} &
&
&
{\small{}Sichuan} &
&
{\small{}Hainan} &
\tabularnewline
&
{\small{}Guizhou} &
&
&
{\small{}Guizhou} &
&
{\small{}Sichuan} &
\tabularnewline
&
{\small{}Shaanxi} &
&
&
{\small{}Yunnan} &
&
{\small{}Gansu} &
\tabularnewline
&
{\small{}Gansu} &
&
&
{\small{}Shaanxi} &
&
{\small{}Qinghai} &
\tabularnewline
&
\textbf{\small{}Yunnan} &
&
&
{\small{}Gansu} &
&
{\small{}Ningxia} &
\tabularnewline
&
{\small{}Ningxia} &
&
&
{\small{}Ningxia} &
&
{\small{}Xinjiang} &
\tabularnewline
&
&
&
&
{\small{}Xinjiang} &
&
&
\tabularnewline
\bottomrule
\end{tabular}
\par\end{centering}
\medskip{}
\raggedright{}\small Note: Provinces in bold highlight different
results from the original paper.
\end{table}
The same as in the original paper, in all three specifications the
information criterion determines two groups. Displayed in Table \ref{tab:chinagdp_result_table},
the classification results by \texttt{Rmosek} are identical to those
by \texttt{CVX} reported in \citet[Table A11]{chinagdp} in the second
and the third specifications. However, in the first specification
where all five indicators are included, we observe deviance across
computing platforms. Liaoning, Zhejiang, Shandong and Yunnan are moved
into group 2 according to \texttt{Rmosek} whereas they are left in
the group 1 in \texttt{CVX} results. Since this is a short $T$ panel
data with $T=8$, the numerical stability is more fragile when we
include more regressors.
A key observation in \citet{chinagdp} is that group 1 gathers Beijing,
Shanghai and Hainan, the three provinces with the highest GDP shares
of the tertiary sector associated with the provision of services.
The \texttt{Rmosek} result retains this feature. As shown in Figure
\ref{fig:Tertiary-industry-GDP}, the tertiary industry GDP shares
of the four provinces excluded from the Beijing-Shanghai-Hainan group
doe not reach the national aggregate share, except for Liaoning in
a single year 2004. Given their relatively lower shares, it is sensible
that \texttt{Rmosek} removes them out of the high-share group.
\begin{figure}
\caption{\label{fig:Tertiary-industry-GDP} GDP Shares of Tertiary Sector}
\medskip{}
\centering{}\includegraphics[width=0.9\textwidth]{teritiary_ratio.png}
\end{figure}
In the parameter tuning process, we compute the information criterion
for $K=1,2,3,4$ and 10 candidate $\lambda$ values.\footnote{We generate $\lambda=c\mathrm{var}\left(y\right)T^{-\frac{1}{3}}$
where the constant $c$ varies from $0.001$ to $0.01$.} We report the CPU time consumed by the parameter tuning process in
Table \ref{tab:Running-Time-china-gdp}. \texttt{Rmosek} is about
$8$ to $9.6$ times faster than \texttt{CVX}. It demonstrates the
speed gain of \texttt{Rmosek} in real applications.
\begin{table}
\caption{\label{tab:Running-Time-china-gdp}Running Time (in second): Replication
of \citet[Table A11]{chinagdp}}
\medskip{}
\begin{centering}
\begin{tabular}{cccc}
\toprule
&
{\small{}With Light} &
{\small{}Without Light} &
{\small{}Without Light and Tax}\tabularnewline
\midrule
\texttt{Rmosek} &
$72.39$ &
$96.37$ &
$40.41$\tabularnewline
\texttt{CVX} &
$692.23$ &
$774.96$ &
$337.10$\tabularnewline
\bottomrule
\end{tabular}
\par\end{centering}
\end{table}
\section{Relaxed Empirical Likelihood \label{sec:Relaxed-Empirical-Likelihood}}
Besides the regression setting in Section \ref{sec:Classifier-Lasso},
convex programming is also useful in structural econometric estimation.
Consider the models with a ``true'' parameter $\beta_{0}$ satisfying
the unconditional moment condition $\mathbb{E}\left[g\left(Z_{i},\beta_{0}\right)\right]=\boldsymbol{0}_{m}$,
where $\left\{ Z_{i}\right\} _{i=1}^{n}$ is the observed data, $\beta\in\mathcal{B\subset}\mathbb{R}^{D}$
is a finite dimensional vector in the parameter space $\mathcal{B}$,
and $g$ is an $\mathbb{R}^{m}$-valued moment function. GMM \citep{hansen1982GMM}
and empirical likelihood (EL) \citep{owen1988empirical,qin1994empirical}
are two workhorses dealing with moment restriction models. In particular,
EL solves
\[
\max_{\beta\in\mathcal{B},\pi\in\Delta_{n}}\,\sum_{i=1}^{n}\log\pi_{i}\quad\text{s.t.}\quad\sum_{i=1}^{n}\pi_{i}g\left(Z_{i},\beta\right)=\boldsymbol{0}_{m}
\]
where $\Delta_{n}=\left\{ \pi\in\left[0,1\right]^{n}:\sum_{i=1}^{n}\pi_{i}=1\right\} $
is the $n$-dimensional probability simplex. However, neither GMM
nor EL can be used to estimate a model with more moment equalities
than observations, i.e. $m>n$. To make the optimization feasible,
\citet{REL} relaxes the equality restriction $\sum_{i=1}^{n}\pi_{i}g_{i}\left(\beta\right)=\boldsymbol{0}_{m}$
in EL. REL is defined as the solution to
\[
\max_{\beta\in\mathcal{B}}\max_{\pi\in\Delta_{n}^{\lambda}\left(\beta\right)}\,\sum_{i=1}^{n}\log\pi_{i}
\]
where
\[
\Delta_{n}^{\lambda}\left(\beta\right)=\left\{ \pi\in\Delta_{n}:\big|\sum_{i=1}^{n}\pi_{i}h_{ij}\left(\beta\right)\big|\leq\lambda,\:j=1,2,\cdots,m\right\}
\]
is a relaxed simplex, $\lambda\geq0$ is a tuning parameter, $h_{ij}\left(\beta\right)=g_{j}\left(Z_{i},\beta\right)/\hat{\text{\ensuremath{\sigma}}}_{j}\left(\beta\right)$,
$g_{j}\left(Z_{i},\beta\right)$ is the $j$-th component of $g\left(Z_{i},\beta\right)$,
and $\hat{\text{\ensuremath{\sigma}}}_{j}\left(\beta\right)$ is the
sample standard deviation of $\left\{ g_{j}\left(Z_{i},\beta\right)\right\} _{i=1}^{n}$.
The formulation of REL is inspired by Dantzig selector \citep{dantzig}.
\bigskip
\begin{example}[Dantzig selector]
Similar to Lasso, Dantzig selector also produces a sparse solution
to the linear regression model. Dantzig selector can be written as
\[
\min_{\beta}\,\lVert\beta\rVert_{1}\ \ \text{s.t.}\ \ \lVert X'\left(y-X\beta\right)\rVert_{\infty}\leq\lambda,
\]
where $\lambda$ is a tuning parameter. We can immediately reformulate
it as a linear programming problem
\begin{gather*}
\min_{\beta^{+},\beta^{-}}\:e'\beta^{+}+e'\beta^{-}\\
\text{s.t.}\ \ X'y-\lambda e\leq\left(X'X\right)\left(\beta^{+}-\beta^{-}\right)\leq X'y+\lambda e\\
\beta^{+},\beta^{-}\geq0.
\end{gather*}
It is readily solvable using the \texttt{R} package \texttt{quantreg}
\citep{quantreg}.
\end{example}
\bigskip
Dantzig selector slacks the sup-norm of the first-order condition
for optimality. REL borrows the idea to estimate a finite-dimensional
parameter in a structural economic model defined by many moment equalities.
Comparing to Dantzig selector, REL uses a nonlinear objective function.
It is still convex (in minus likelihood) but \texttt{quantreg} that
deals with linear programming problems is no longer applicable.
Similar to standard EL, REL's optimization involves an inner loop
and an outer loop. The outer loop for $\beta$ is a general low-dimensional
nonlinear optimization, which can be solved by Newton-type methods.
With the linear constraints and the logarithm objective, the inner
loop is convex in $\pi=\left(\pi_{i}\right)_{i=1}^{n}$. For each
$\beta$, the inner problem can be formulated as a \emph{separable
convex optimization problem} in the matrix form
\begin{gather*}
\max_{\pi}\,\sum_{i=1}^{n}\log\pi_{i}\\
\text{s.t.}\ \ \begin{bmatrix}1\\
-\lambda\\
\vdots\\
-\lambda
\end{bmatrix}\leq\begin{bmatrix}1 & 1 & \cdots & 1\\
h_{11}\left(\beta\right) & h_{21}\left(\beta\right) & \cdots & h_{n1}\left(\beta\right)\\
\vdots & \vdots & \ddots & \vdots\\
h_{1m}\left(\beta\right) & h_{2m}\left(\beta\right) & \cdots & h_{nm}\left(\beta\right)
\end{bmatrix}\begin{bmatrix}\pi_{1}\\
\pi_{2}\\
\vdots\\
\pi_{n}
\end{bmatrix}\leq\begin{bmatrix}1\\
\lambda\\
\vdots\\
\lambda
\end{bmatrix}\\
0\leq\pi_{i}\leq1,\:\text{for each}\,i=1,2,\cdots,n
\end{gather*}
and it is readily solvable in \texttt{Rmosek} by translating the mathematical
expression into computer code.
\subsection{Replication}
We follow the simulation design in \citet[Section 4]{REL}, which
is described in Appendix \ref{subsec:REL-DGP}. Table \ref{tab:relsimulation}
reports the bias and RMSE of the estimation of $\hat{\beta}_{1}$,
implemented purely in \texttt{R} with the inner loop by \texttt{Rmosek}
and the outer loop by \texttt{nloptr}. The results are close to those
in \citet{REL}, where the code is written in \texttt{MATLAB} with
the outer loop handled by the function \texttt{fmincon} and the inner
loop by \texttt{CVX} solved by \texttt{MOSEK}.
\begin{table}[h]
\caption{\label{tab:relsimulation}Estimation of $\beta_{1}$ in linear IV
model with REL: Replication of \citet{REL} }
\medskip{}
\centering{}
\begin{tabular}{rcccccc}
\toprule
&
&
\multicolumn{2}{c}{Replication} &
&
\multicolumn{2}{c}{Original Results}\tabularnewline
\cmidrule{3-4} \cmidrule{4-4} \cmidrule{6-7} \cmidrule{7-7}
$\left(n,m\right)$ &
&
Bias &
RMSE &
&
Bias &
RMSE\tabularnewline
\midrule
$\left(120,80\right)$ &
&
-0.020 &
0.135 &
&
-0.004 &
0.113\tabularnewline
$\left(120,160\right)$ &
&
-0.018 &
0.162 &
&
-0.012 &
0.143\tabularnewline
$\left(240,80\right)$ &
&
-0.004 &
0.078 &
&
-0.006 &
0.071\tabularnewline
$\left(240,160\right)$ &
&
-0.008 &
0.093 &
&
-0.009 &
0.077\tabularnewline
\bottomrule
\end{tabular}
\end{table}
We also experiment with other numerical alternatives. Since the scale
of the optimization problems here is much smaller than C-Lasso, the
inner loop can be correctly solved by \texttt{Rmosek}, \texttt{CVXR},
\texttt{CVX} in \texttt{MATLAB}, or even \texttt{nloptr}. These four
methods produce virtually identical inner loop results up to rounding
errors. This finding confirms the robustness of the \texttt{R} environment
in high-dimensional optimization. The difference in Table \ref{tab:relsimulation},
therefore, is attributed to the outer loop between the function \texttt{nloptr}
in \texttt{R} and the function \texttt{fmincon} in \texttt{MATLAB}.
\begin{table}
\caption{\label{tab:rel_running_time}Running time of REL's inner loop \textcolor{black}{(in
second)}}
\medskip{}
\centering{}
\begin{tabular}{crrrr}
\toprule
$\left(n,m\right)$ &
$\left(120,80\right)$ &
$\left(120,160\right)$ &
$\left(240,80\right)$ &
$\left(240,160\right)$\tabularnewline
\midrule
\texttt{Rmosek} &
2.995 &
4.378 &
10.510 &
17.206\tabularnewline
\texttt{nloptr} &
64.904 &
117.533 &
115.738 &
226.661\tabularnewline
\texttt{CVXR} &
31.241 &
43.909 &
42.435 &
136.007\tabularnewline
\texttt{CVX} &
41.441 &
54.095 &
65.846 &
88.982\tabularnewline
\bottomrule
\end{tabular}
\end{table}
To evaluate the computational cost, we record the time spent in the
inner loop. With $100$ sets of data generated by the same DGP for
each sample size, we fix $\beta=\left(0.9,\,0.9\right)$ and only
numerically solve the inner loop. Since four approaches have virtually
identical inner loop results, we only report the running time of each
method in Table \ref{tab:rel_running_time}. Although \texttt{CVXR}
and \texttt{nloptr} are able to correctly solve the problem thanks
to its small scale, \texttt{Rmosek} remains $4$ to $30$ times faster
than these alternatives. We conjecture that bigger speed gain would
be observed in a problem of larger scale.
\section{Conclusion\label{sec:Conclusion}}
In this note, we demonstrate numerical implementation via \texttt{Rmosek}
of two examples of high-dimensional econometric estimators. The convenience
and reliability of high-dimensional convex optimization in \texttt{R}
will open new possibilities to create estimation procedures. In the
era of big data, we are looking forward to witnessing more algorithms
blossoming and flourishing along with theoretical research of high-dimensional
models.
\bigskip
\singlespacing\bibliographystyle{chicagoa}
\bibliography{sample}
\onehalfspacing
\newpage{}