% verifier version packages
% coef(summary) pb philosophique
% refaire phrase sur Vuong dans pscl
% DF pour le summary
% Revoir le logLik (AIC)
% SMITH:03 -> SMITH:02

\documentclass[nojss]{jss}
\usepackage{rotating}
\usepackage{amsmath}
\usepackage{wasysym}
%\usepackage{threeparttable}
\usepackage{amssymb, amsfonts}
\usepackage{longtable}
\usepackage{lscape}
\newcommand{\change}[1]{#1}

%\VignetteIndexEntry{Multiple Hurdle Models in R: the mhurdle Package}
%\VignetteDepends{Formula, truncreg, maxLik}
%\VignetteKeywords{Limited dependent variable, tobit, hurdle models, econometric computing, R}
%\VignettePackage{mhurdle}

\title{Multiple Hurdle Tobit Models in \proglang{R}: The \pkg{mhurdle} Package}

\Plaintitle{Multiple Hurdle Tobit Models in R: The mhurdle Package}

% \footnote{This package has been developed as part
% of a PhD dissertation carried out by St\'ephane Hoareau
% \cite{Hoareau/09} at the University of La R\'eunion under the
% supervision of Fabrizio Carlevaro and Yves Croissant.}


\author{Fabrizio Carlevaro\\Universit\'e de Gen\`eve \And
 Yves Croissant\\Universit\'e de la R\'eunion \And
 St\'ephane Hoareau\\Universit\'e de la R\'eunion}

\Plainauthor{Fabrizio Carlevaro, Yves Croissant, St\'ephane Hoareau}

\Address{
Fabrizio Carlevaro\\
Facult\'e des sciences \'economiques et sociales\\
Universit\'e de Gen\`eve\\
Uni Mail\\
40 Bd du Pont d'Arve\\
CH-1211 Gen\`eve 4\\
Telephone: +41/22/3798914\\
E-mail:\email{fabrizio.carlevaro@unige.ch}
\\
\\
Yves Croissant\\
Facult\'e de Droit et d'Economie\\
Universit\'e de la R\'eunion\\
15, avenue Ren\'e Cassin\\
BP 7151\\
F-97715 Saint-Denis Messag Cedex 9\\
Telephone: +33/262/938446\\
E-mail: \email{yves.croissant@univ-reunion.fr}
\\
\\
St\'ephane Hoareau\\
Facult\'e de Droit et d'Economie\\
Universit\'e de la R\'eunion\\
15, avenue Ren\'e Cassin\\
BP 7151\\
F-97715 Saint-Denis Messag Cedex 9\\
Telephone: +33/262/938446\\
E-mail: \email{stephane.hoareau@univ-reunion.fr}
}

%% need no \usepackage{Sweave.sty}

\Abstract{ \pkg{mhurdle} is a package for \proglang{R} enabling the
  estimation of a wide set of regression models where the dependent
  variable is left censored at zero, which is typically the case in
  household expenditure surveys.  These models are of particular
  interest to explain the presence of a large proportion of zero
  observations for the dependent variable by means of up to three
  censoring mechanisms, called hurdles. For the analysis of censored
  household expenditure data, these hurdles express a good selection
  mechanism, a desired consumption mechanism and a purchasing
  mechanism, respectively. \pkg{mhurdle} models are specified in a
  fully parametric form and estimated using the maximum likelihood
  method for random samples. Model evaluation and selection are
  tackled by means of goodness of fit measures and Vuong
  tests. Software rationale and user's guidelines are presented and
  illustrated with a real-world example.}

\Keywords{Households' expenditure survey analysis, censored regression
  models, hurdle models, Tobit models, maximum likelihood estimation,
  goodness of fit measures, Vuong tests, \proglang{R}}

\Plainkeywords{Households' expenditure survey analysis, censored
  regression models, hurdle models, Tobit models, maximum likelihood estimation,
  goodness of fit measures, Vuong tests, R}

\begin{document}


\maketitle




\section{Introduction}
<<echo=FALSE,results=hide, label = init>>=
options(prompt= "R> ", useFancyQuotes = FALSE)
@
Data collected by means of households' expenditure survey may present
a large proportion of zero expenditures due to many households
recording, for one reason or another, no expenditure for some
items. Analyzing these data requires to model any expenditure with a
large proportion of nil observations as a dependent variable left
censored at zero.

Since the seminal paper of \citet{TOBIN/58}, a large econometric
literature has been developed to deal correctly with this problem of
zero observations. The problem of censored data has also been treated
for a long time in the statistics literature dealing with survival
models.

In applied microeconometrics, different decision mechanisms have been
put forward to explain the appearance of zero expenditure
observations. The original Tobin's model takes only one of these
mechanisms into account. With \pkg{mhurdle}, up to three mechanisms
generating zero expenditure observations may be introduced in the
model \footnote{This package is an improved version of a package first
  developed as part of a PhD dissertation carried out by St\'ephane
  \citet{HOAR/09} at the University of R\'eunion under the supervision
  of Fabrizio Carlevaro and Yves Croissant.}. More specifically, we
consider the following three zero expenditure generating mechanisms.

\begin{description}
\item[A good selection mechanism (hurdle 1)]. According to this
  mechanism, the consumer\footnote{The consumer we are reffering to is
    that of the microeconomic theory, an abstract economic agent
    responsible of the decisions of a consumption unit that may be an
    individual, a family, a household. According to the economic
    litterature, we term this concept ``the consumer'' by
    convenience.} first decides which goods to include in its choice
  set and, as a consequence, he can discard some marketed goods
  because he dislikes them (like meat for vegetarians or wine for
  non-drinkers) or considers them harmful (like alcohol, cigarettes,
  inorganic food, holidays in dangerous countries), among
  others.\\
  This censoring mechanism has been introduced in empirical demand
  analysis by \citet{CRAGG/71}. It allows to account for the
  non-consumption of a good as a consequence of a fundamentally
  non-economic decision motivated by ethical, psychological or social
  considerations altering the consumer's preferences.

\item[A desired consumption mechanism (hurdle 2)]. According to this
  mechanism, once a good has been selected, the consumer decides which
  amount to consume and, as a consequence of his preferences,
  resources and selected good prices, its rational decision can turn
  out to be a
  negative desired consumption level leading to a nil consumption.\\
  The use of this mechanism, to explain the presence
  of zero observations in family expenditure surveys, was introduced
  by \citet{TOBIN/58}. Its theoretical relevance has been later
  rationalised by the existence of corner solutions to the
  microeconomic problem of rational choice of the neoclassical
  consumer. See section 10.2 of \citet{AMEM/85}, for an elementary
  presentation of this issue, and chapter 4 of \citet{PUDNEY/89}, for
  a more comprehensive one.

\item[A purchasing mechanism (hurdle 3)]. According to this mechanism,
  once a consumption decision has been taken, the consumer sets up the
  schedule at which to buy the good and, as a consequence of its
  purchasing strategy, zero expenditure may be observed if the survey
  by which these data are collected is carried out over a too short
  period with respect to the frequency at which the good is bought.\\
  This censoring mechanism has been introduced in empirical demand
  analysis by \citet{DEATO/IRISH/84}. It allows to account for the
  non-purchase of a good not because the good is not consumed but
  because it is a durable or a storable good infrequently bought. By
  the same token, this mechanism allows to derive from observed
  expenditures, the rate of use of a durable good or the rate of
  consumption of a stored non durable good.
\end{description}

For each of these censoring mechanisms, a continuous latent variable
is defined, indicating that censoring is in effect when the latent
variable is negative. These latent variables are modelled as the sum of
a linear function of explanatory variables and of a normal random
disturbance, with a possible correlation between the disturbances of
different latent variables in order to account for a 
possible simultaneity of the decisions modelled by censoring mechanisms.
\change{To model possible departures of the observed dependent 
variable to normality, we use flexible transformations allowing 
to rescale skewed or leptokurtic random variables to normality.}   
By combining part or the whole set of
these censoring mechanisms, we generate a set of non-nested parametric
models that can be used to explain censored expenditure data depending
on the structural censoring mechanisms that a priori information
suggests to be at work.

These formal models have been primarily developed to deal with
censored household expenditure data, and numerous applications have
been carried out in this field. Complementing a previous survey by
\citet{SMITH:02}, Table~\ref{tab:exconso} gives an updated overview of these
studies. We note a late popularity of Cragg's approach, 
as the first applications of hurdle models are published in the late of 
1980s, namely almost two decades after the publication of Cragg's paper
in Econometrica. Since then, a large variety of demand models including 
one or two among the previous three censoring mechanisms are estimated.
However, none of these studies use the three censoring mechanisms 
we consider, jointly. From the 1990s on, many studies account for deviations 
of the desired consumption relation to homoscedasticity and normality by
modelling the standard error of this variable as a non negative 
parametric function of some explanatory variables and by transforming 
its distribution to normality using either the Box-Cox or the inverse
hyperbolic sine transformations. The estimation of a correlation coefficient 
between disturbances is also performed in several of these studies, with 
an increasing succes over time, in terms of statistical significance of estimates.

\input{biblio.tex}
% revoir les \cite ci-dessous et le WODJ:2020

The practical scope of multiple hurdle Tobit models is not restricted
to empirical demand analysis but has been fruitfully used in other
fields of economics. This includes labor
economics~\citep[][]{ELEK:KOLL:REIZ:SZAB:11}, contingent valuation
studies \citep[][]{SAZS:RAUS:08,MART:06}, finance
\citep[][]{MOFF:05}, sport activities \citep[][]{HUMP:RUSE:10},
internet use \citep[][]{WODJ:2020}, gambling
\citep[][]{HUMP:LEE:SOEB:09}, production \citep[][]{AKPA:NKAN:ESSI:12,
  MAL:ANIK:BAUE:SCHM:12, OKEL:KIRU:GITO:12, TEKL:DADI:YAMI:DANA:06}

%% It is worth mentioning that, although this formal model has been
%% primarily developed to deal with censored household expenditure data,
%% its practical scope is not restricted to empirical demand analysis. A
%% quite natural other area of application is represented by the
%% empirical analysis of labour supply. In this context, hurdle 1 can
%% indeed be reinterpreted as a non-economic mechanism of labour market
%% participation; hurdle 2 as a desired working hours mechanism based on
%% the neoclassical model of labour supply that can generate negative
%% desired working hours leading to a nil labour supply; hurdle 3 as an
%% unemployment mechanism explaining zero hours worked as a result of
%% spells of unemployment. Note also that even within the realm of demand
%% analysis, the economic interpretation of hurdles 1, 2 and 3 may
%% require to be adapted to the specific features of available data, as
%% we illustrate by an empirical application presented at the end of the
%% paper (see section \ref{sec:examples}).

Our hurdle models are specified as fully parametric models allowing
estimation and inference within an efficient maximum likelihood
framework. In order to identify a relevant model specification,
goodness of fit measures for model evaluation and selection, as well
as Vuong tests for discriminating between nested, strictly non nested
and overlapping models have been implemented in \pkg{mhurdle}
package. Vuong tests remarkably permit to compare two competing models
when both, only one, or neither of them contain the true mechanism
generating the sample of observations.  More precisely, such tests
allow to assess which of the two competing models is closest to the
true unknown model according to the Kullback-Leibler information
criterion. Therefore, such symmetric tests are not intended, as
classical Neyman-Pearson tests, to pinpoint the chimeric true model,
but to identify a best parametric model specification (with respect to
available observations) among a set of competing specifications. As a
consequence, they can provide inconclusive results, which prevent from
disentangling some competing models, and when they are conclusive,
they don't guarantee an identification of the relevant model
specification.

Survival models are implemented in \proglang{R} with the
\pkg{survival} package of \citet{SURVA/08}. It has also close links
with the problem of selection bias, for which some methods are
implemented in the \pkg{sampleSelection} package of
\citet{TOOME/HENNI/08}. It is also worth mentioning that a convenient
interface to \code{survreg}, called \code{tobit}, particularly aimed
at econometric applications is available in the \pkg{AER} package of
\citet{KLEI/ZEIL/08}. More enhanced censored regression models (left
and right censoring, random effect models) are available in the
\pkg{censReg} package \citep[][]{HENN:13}. Some flavor of hurdle
models have also been developed for count data and are implemented by
the \code{hurdle} of the \pkg{pscl} package
\citep[][]{ZEIL:KLEI:JACK:08}.

The paper is organised as follows: Section \ref{sec:limdepvar}
presents the rationale of our modelling strategy.  Section
\ref{sec:estevalsel} presents the theoretical framework for model
estimation, evaluation and selection. Section \ref{sec:software}
discusses the software rationale used in the package. Section
\ref{sec:examples} illustrates the use of \pkg{mhurdle} with a
real-world example. Section \ref{sec:conclusion} concludes.

\section{Modelling strategy}
\label{sec:limdepvar}

\subsection{Model specification}
\label{sec:modelspec}

Our modelling strategy is intended to model the level $y$ of expenditures
of a household for a given good or service during a given period of observation.
To this purpose, we use up to three zero expenditure generating mechanisms,
called hurdles, and a demand function.

Each hurdle is represented by a probit model resting on one of the following three
latent dependent variables relations:

\begin{equation}\label{eq1}
\left\{
\begin{array}{l}
  y_1^* = \beta_1^\top x_1 + \epsilon_1 \\[4pt]
  y_2^* = \beta_2^\top x_2 + \epsilon_2 \\[4pt]
  y_3^* = \beta_3^\top x_3 + \epsilon_3 \\
\end{array}
\right.
\end{equation}

where $x_1$, $x_2$, $x_3$ stand for column-vectors of explanatory
variables (called covariates in the followings), $\beta_1$, $\beta_2$,
$\beta_3$ for column-vectors of the impact coefficients of the
explanatory variables on the continuous latent dependent variables $y_1^*$, $y_2^*$,
$y_3^*$ and $\epsilon_1$, $\epsilon_2$, $\epsilon_3$ for normal random
disturbances. Since variables $y_1^*$ and $y_3^*$ are never observed, 
contrary to $y_2^*$, the units of measurement of $\epsilon_1$ and $\epsilon_3$ are
not identified. Hence, these disturbances are normalized by setting their variances 
equal to 1, i.e. by identifying them to standard normal random variables.

\begin{itemize}

\item Hurdle 1 models the household decision of selecting or not
  selecting the good we consider as a relevant consumption good,
  complying with household's ethical, psychological and social
  convictions and habits. This good selection mechanism explains the
  outcome of a binary choice that can be coded by a binary variable
  $I_1$ taking value 1 if the household decides to enter the good in
  its basket of relevant consumption goods and 0 otherwise. The
  outcome of this binary choice is modelled by associating the decision
  to select the good to positive values of the latent variable $y_1^*$
  and that to reject the good to negative values of $y_1^*$.
  Therefore, good selection or rejection is modelled as a probability
  choice where selection occurs with probability $P(I_1=1)=P(y_1^*>0)$
  and rejection
  with probability $P(I_1=0)=P(y_1^*\leq 0)=1-P(y_1^*>0)$.\\
  Note that if this mechanism is inoperative, this probit model must
  be replaced by a singular probability choice model where
  $P(I_1=1)=1$ and $P(I_1=0)=0$.

\item Hurdle 2 models the household decision of consuming or not
  consuming the selected good, given its actual economic
  conditions. This desired consumption mechanism explains the outcome
  of a binary choice coded by a binary variable $I_2$ taking value 1
  if the household decides to consume the good and 0 otherwise.  The
  outcome of this binary choice is modelled by associating the
  decision to consume the selected good to a positive value of its
  desired consumption level, represented by the latent variable
  $y_2^*$, and that of not to consume the good to negative values of
  $y_2^*$. Therefore, when this zero expenditure generating mechanism
  is operative, it also models the level of desired consumption
  expenditures by means of a Tobit model identifying the desired
  consumption expenditures to the value of latent variable $y_2^*$,
  when it is positive, and to zero, when it is negative.\\
  Conversely, when the desired consumption mechanism is inoperative,
  implying that the desired consumption cannot be a corner solution of
  a budget constrained problem of utility minimisation, we must
  replace not only the probit model explaining the variable $I_2$ by a
  singular probability choice model where $P(I_2=1)=1$, but also the
  Tobit demand function by a demand model enforcing non-negative
  values on the latent variable $y_2^*$.  \citet{CRAGG/71}
  suggested two types of functional forms for this demand model, namely
  a log-normal functional form :

  \begin{equation}\label{eq2}
  \ln y_2^* = \beta_2^\top x_2+\epsilon_2
  \end{equation}

  and a truncated normal functional form where $y_2^*$ is generated by
  a linear relationship $y_2^*=\beta_2^\top x_2 + \epsilon_2$ with
  $\epsilon_2$ distributed according to a normal distribution
  left-truncated at $\epsilon_2 = -\beta_2^\top x_2$.  Nevertheless,
  to avoid a cumbersome analytic presentation of our models, in the
  following we only consider the log-normal model specification.
  More flexible generalizations of these functional
  forms will be discussed in section 2.3.

\item Hurdle 3 models the household decision to purchase or not to
  purchase the good during the survey period over which expenditure
  data are collected. This purchasing mechanism also explains the
  outcome of a binary choice, coded by a binary variable $I_3$ taking
  value 1 if the household decides to buy the good during the period
  of statistical observation and 0 otherwise.  The probit model we use
  associates the purchasing decision to positive values of latent
  variable $y_3^*$ and that of not purchasing to negative
  values of $y_3^*$.\\
  By assuming that consumption and purchases are uniformly distributed
  over time, but according to different timetables entailing a
  frequency of consumption higher than that of purchasing, we can also
  interpret the probability $P(I_3=1)=P(y_3^*>0)$ as measuring the
  share of purchasing frequency to that of consumption during the
  observation period. This allows to relate the observed level of
  expenditures $y$ to the unobserved level of consumption $y_2^*$
  during the observation period, using the following identity:

  \begin{equation}\label{eq3}
  y = \frac{y_2^*}{P(I_3=1)} I_1 I_2 I_3 .
  \end{equation}

  When the purchasing mechanism is inoperative, the previous probit
  model must be replaced by a singular probability choice model where
  $P(I_3=1)=1$.  In such a case, the observed level of expenditures is
  identified to the level of consumption, implying $y=y_2^* I_1 I_2$.

\end{itemize}

A priori information (theoretical or real-world knowledge) may suggest
that one or more of these censoring mechanisms are not in effect. For
instance, we know in advance that all households purchase food
regularly, implying that the first two censoring mechanisms are
inoperative for food.  In this case, the relevant model is defined by
only two relations: one defining the desired consumption level of
food, according to a log-normal or a truncated normal specification,
and the other the decision to purchase food during the observation
period.

Figure \ref{fig:arbre2} outlines the full set of special models that
can be generated by selecting which of these three mechanisms are in
effect and which are not. It shows that 8 different models can be
dealt with by means of \pkg{mhurdle} package.  To use a mnemonically
rule, we number the models by 3 binary digits, each of which indicates
if a censoring mechanism is or is not in effect, using figures
\code{1} and \code{0}, respectively.  For example, \code{011}
indicates the model for which hurdle 1 (good selection mechanism) is
not in effect while hurdles 2 (desired consumption mechanism) and 3
(purchasing mechanism) are.


\begin{sidewaysfigure}
\caption{\label{fig:arbre2} The full set of mhurdle special models.}
  \includegraphics[width=20cm,height=17cm]{typo1.pdf}
\end{sidewaysfigure}

Among these models, one is not concerned by censored data, namely
model \code{000}. This model is relevant only for modelling uncensored
samples. All the other models are potentially able to analyse censored
samples by combining up to the three censoring mechanisms described
above. With the notable exception of the standard Tobit model \code{010}, that
can be estimated also by the \pkg{survival} package of
\citet{SURVA/08} or the \pkg{AER} package of \citet{KLEI/ZEIL/08}, these
models cannot be found in an other \proglang{R} package.

Some of \pkg{mhurdle} models have already been used in applied
econometric literature. In particular, model \code{100} is a
single-hurdle good selection model originated by \citet{CRAGG/71}
by assuming independence between disturbances $\epsilon_1$ and $\epsilon_2$.
\change{The dependent version of this model
may be viewed as a sample selection model in which only the desired consumption 
is observed, but it diffears from this model
populazied by Heckman (1979) to illustrate linear regression model estimation
given sample selection, in that desired consumption is generated by a demand model
enforcing non negative values on latent variable $y_2^*$. 
Still, in many applications presented in Table~\ref{tab:exconso},
Heckman's sample selection model is used as a dependent single-hurdle
good selection model in which good selection decision is assumed to dominate good
consumption decision. From our point of view, this model is theoretically
misspecified to analyse a latent dependent variable that could take negative values 
while assuming that hurdle 2 is not in effect.} 

The double-hurdle model 
\code{110} combining independent good selection (hurdle 1) and desired
consumption (hurdle 2) censoring mechanisms is also due to
\citet{CRAGG/71}. An extension of this double-hurdle model to
dependent censoring mechanisms has been originated by
\citet{BLUNDELL/87}.

P-Tobit model \code{011} is due to \citet{DEATO/IRISH/84} and explains zero
purchases by combining the desired consumption censoring mechanism
(hurdle 2) with the purchasing censoring mechanism (hurdle 3).
Model \code{001} is a single-hurdle model not yet used in applied demand analysis,
where the censoring mechanism in effect is that of infrequent purchases
(hurdle 3).

Among the original models encompassed by \pkg{mhurdle}, models \code{101} is
a double-hurdle model combining good selection (hurdle 1) and
purchasing (hurdle 3) mechanisms to explain censored samples. Model
\code{111} is an original triple-hurdle model originated in \citet{HOAR/09}.
This model explains censored purchases either as the result of good
rejection (hurdle 1), negative desired consumption (hurdle 2) or
infrequent purchases (hurdle 3).

To derive the form of the probability distribution of the observable
dependent variable $y$, we must specify the joint distribution of the
random disturbances entering the structural relations of these models.

\begin{itemize}

\item Models \code{111} and \code{101} are trivariate hurdle models as they involve
  disturbances $\epsilon_1$, $\epsilon_2$ and $\epsilon_3$,
  distributed according to the trivariate normal density function:

  \begin{equation}\label{eq4}
    \frac{1}{\sigma}\phi\left(\epsilon_1,\frac{\epsilon_2}{\sigma},\epsilon_3;
      \rho_{12},\rho_{13},\rho_{23}\right),
  \end{equation}

  where

  $$
  \phi(z_1,z_2,z_3;\rho_{12},\rho_{13},\rho_{23})={\displaystyle
    \frac{\exp\left\{-\frac{\rho^{11}z_1^2+\rho^{22}z_2^2+\rho^{33}z_3^2
          -2[\rho^{12}z_1z_2+\rho^{13}z_1z_3+\rho^{23}z_2z_3]}{2}\right\}}
    {\sqrt{(2\pi)^3}\mid R\mid}},
  $$
  with
  $$
  |R|=1-\rho_{12}^2-\rho_{13}^2-\rho_{23}^2+2\rho_{12}\rho_{13}\rho_{23},
  $$
  $$
  \rho^{11}=\frac{1-\rho_{23}^2}{\mid R\mid},\quad
  \rho^{22}=\frac{1-\rho_{13}^2}{\mid R\mid},\quad
  \rho^{33}=\frac{1-\rho_{12}^2}{\mid R\mid},
  $$
  $$
  \rho^{12}=\frac{\rho_{12}-\rho_{13}\rho_{23}}{\mid R\mid},\quad
  \rho^{13}=\frac{\rho_{13}-\rho_{12}\rho_{23}}{\mid R\mid},\quad
  \rho^{23}=\frac{(\rho_{23}-\rho_{12}\rho_{13})}{\mid R\mid},
  $$

  denotes the density function of a standard trivariate normal
  distribution and $\rho_{12}$, $\rho_{13}$, $\rho_{23}$ the
  correlation coefficients between the couples of normal standard
  random variables $z_1$ and $z_2$, $z_1$ and $z_3$, $z_2$ and $z_3$,
  respectively.

\item Models \code{011} and \code{001} are bivariate hurdle models as they involve disturbances
  $\epsilon_2$ and $\epsilon_3$, distributed according to the bivariate normal
  density function:

  \begin{equation}\label{eq5}
  \frac{1}{\sigma}\phi\left(\frac{\epsilon_2}{\sigma},\epsilon_3;\rho_{23}\right),
  \end{equation}

  where

  $$
  \phi(z_1,z_2;\rho)=\frac{\exp\left\{-\frac{z_1^2+z_2^2-2\rho z_1z_2}{2(1-\rho^2)}\right\}}
  {2\pi\sqrt{1-\rho^2}}
  $$

  denotes the density function of a standard bivariate normal distribution with
  correlation coefficient $\rho$.

\item Models \code{110} and \code{100} are also bivariate hurdle models but they involve disturbances
  $\epsilon_1$ and $\epsilon_2$ which density function is therefore written as:

  \begin{equation}\label{eq6}
  \frac{1}{\sigma}\phi\left(\epsilon_1,\frac{\epsilon_2}{\sigma};\rho_{12}\right).
  \end{equation}

\item Finally, models \code{010} and \code{000} are univariate hurdle models involving only
  disturbance $\epsilon_2$, which density function writes therefore:

  \begin{equation}\label{eq7}
  \frac{1}{\sigma}\phi\left(\frac{\epsilon_2}{\sigma}\right),
  \end{equation}

  where

  $$
  \phi(z_1)=\frac{\exp\left\{-\frac{z_1^2}{2}\right\}}{\sqrt{2\pi}}
  $$

  denotes the density function of a standard univariate normal
  distribution.

\end{itemize}

While the assumption of correlated disturbances is intended to account
for the interdependence between latent variables $y_1^*$, $y_2^*$ and
$y_3^*$ unexplained by covariates $x_1$, $x_2$ and $x_3$, a priori
information (theoretical or real-world knowledge) may also suggest to
set to zero some or all correlations between the random disturbances
entering these models, entailing a partial or total independence
between model relations. The use of this a priori information
generates, for each trivariate or bivariate hurdle model of Figure 1,
a subset of special models all nested within the general model from
which they are derived. For a trivariate hurdle model the number of
special models so derived is equal to 7, but for a bivariate hurdle
model only one special model is generated, namely the model obtained
by assuming the independence between the two random disturbances of
the model.

In the following, we shall work out the distribution of our hurdle
models in their general case, but considering the difficulties of
implementing trivariate hurdle models in their full generality, for
these models only the special cases of independence or dependence
between one of hurdles 1 or 3 and the desired consumption equation,
which seems the most relevant for empirical applications, have been
programmed in the present version of \pkg{mhurdle}.  To identify the
presence or absence of assumed dependence between couples of
disturbances of a given \pkg{mhurdle} special model, we add to the 3
binary digit number of the model a letter \code{i} , if independence
is assumed, and a letter \code{d} , otherwise.  For exemple,
\code{101dii} indicates a trivariate model \code{101} for which the
couple of disturbances $(\epsilon_1 ,\epsilon_2)$ are assumed to be
correlated, while $(\epsilon_1 ,\epsilon_3)$ and $(\epsilon_2
,\epsilon_3)$ are not.

\subsection{Likelihood function}

As for the standard Tobit model, the probability distribution of the
observed censored variable $y$ of our hurdle models is a
discrete-continuous mixture, which assigns a probability mass $P(y=0)$
to $y=0$ and a density function $f_+(y)$ to any $y>0$, with:

\begin{equation}\label{eq8}
P(y=0)+\int_0^\infty f_+(y)dy=1.
\end{equation}

The probability mass $P(y=0)=1-P(y>0)$ may be computed by integrating
the joint density function of the latent variables entering the hurdle
model over their positive values.

\begin{itemize}

\item For trivariate hurdle model \code{111}, using the change of
  variables:

  \begin{equation}\label{eq9}
  \left\{
  \begin{array}{l}
   z_1  =  y_1^* - \beta_1^\top x_1 \\[4pt]
   z_2  =  {\displaystyle\frac{y_2^* - \beta_2^\top x_2}{\sigma}} \\[8pt]
   z_3  =  y_3^* - \beta_3^\top x_3 \\
  \end{array}
  \right.
  \end{equation}

  this approach leads to:

  \begin{equation}\label{eq10}
  \begin{array}{l}
  P(y=0) ={\displaystyle 1-\int_{-\beta_1^\top x_1}^\infty \int_{-\frac{\beta_2^\top x_2}{\sigma}}^\infty
  \int_{-\beta_3^\top x_3}^\infty \phi(z_1,z_2,z_3;\rho_{12},\rho_{13},\rho_{23})
  dz_1 dz_2 dz_3}\\[16pt]
  \phantom{P(y=0)} ={\displaystyle 1-\Phi(\beta_1^\top x_1,\frac{\beta_2^\top x_2}{\sigma},
  \beta_3^\top x_3;\rho_{12},\rho_{13},\rho_{23})},
  \end{array}
  \end{equation}

  where $\Phi(z_1,z_2,z_3;\rho_{12},\rho_{13},\rho_{23})$ denotes the distribution
  function of a standard trivariate normal distribution with correlation coefficients
  $\rho_{12}$, $\rho_{13}$ and $\rho_{23}$.

\item For trivariate hurdle model \code{101}, using the change of
  variables:

  \begin{equation}\label{eq11}
  \left\{
  \begin{array}{l}
   z_1 = y_1^* - \beta_1^\top x_1 \\[4pt]
   z_2 = {\displaystyle\frac{\ln y_2^* - \beta_2^\top x_2}{\sigma}} \\[8pt]
   z_3 = y_3^* - \beta_3^\top x_3 \\
  \end{array}
  \right.
  \end{equation}

  this approach leads to:

  \begin{equation}\label{eq12}
  \begin{array}{l}
  P(y=0)={\displaystyle 1-\int_{-\beta_1^\top x_1}^{\infty} \int_{-\infty}^\infty
  \int_{-\beta_3^\top x_3}^{\infty} \phi(z_1,z_2,z_3;\rho_{12},\rho_{13},\rho_{23})
  dz_1 dz_2 dz_3}\\[16pt]
  \phantom{P(y=0)}={\displaystyle 1-\Phi(\beta_1^\top x_1,\beta_3^\top x_3;\rho_{13})},
  \end{array}
  \end{equation}

  where $\Phi(z_1,z_2;\rho)$ denotes the distribution function of a standard
  bivariate normal distribution with correlation coefficient $\rho$.

\item The probability mass $P(y=0)$ for bivariate hurdle models
  \code{011} and \code{110} and univariate hurdle model \code{010} can
  be derived from that of trivariate model \code{111} by eliminating
  hurdles 1, 3, 1 and 3, respectively.  Likewise, this probability for
  bivariate hurdle models \code{100} and \code{001} can be derived
  from that of trivariate hurdle model \code{101} by eliminating
  hurdles 1 and 3, respectively.  Corresponding formulas of $P(y=0)$
  for all this special cases implemented in \proglang{R} are presented
  in Table \ref{tab:typo}, using the following notations:
  $$
  \Phi_1=\Phi(\beta_1^\top x_1),\quad \Phi_2=\Phi\left(\frac{\beta_2^\top x_2}{\sigma}\right),
  \quad \Phi_3=\Phi(\beta_3^\top x_3),
  $$
  $$
  \Phi_{12}=\left(\beta_1^\top x_1,\frac{\beta_2^\top x_2}{\sigma};\rho_{12}\right),\quad
  \Phi_{23}=\left(\frac{\beta_2^\top x_2}{\sigma},\beta_3^\top x_3;\rho_{23}\right),
  $$
  where $\Phi(z)$ denotes the distribution function of a standard
  univariate normal distribution.

\end{itemize}

\begin{sidewaystable}
  \caption{Characteristics of mhurdle special models implemented in
  \proglang{R}\label{tab:typo}}\vspace{.5cm}

  \begin{tabular}{|l|lll|lll|l|l|l|}\hline\hline
id & $h_1$ & $h_2$ & $h_3$ & $\rho_{12}$ & $\rho_{13}$ & $\rho_{23}$  & $P(y=0)$ & $f_+(y)$ & $\mbox{E}(y \mid y > 0)$\\\hline\hline
000 &$\square$&$\square$&$\square$&&&& $0$ & $\frac{1}{\sigma y}\phi\left(\frac{\ln y - \beta_2^{\top}x_2}{\sigma}\right)$ & $\exp\left\{\beta_2^{\top}x_2 + \frac{\sigma^2}{2}\right\}$\\
100i &$\blacksquare$&$\square$&$\square$&$\square$&&& $1 - \Phi_1$& $\frac{1}{\sigma y}\phi\left(\frac{\ln y -\beta_2^{\top}x_2}{\sigma}\right)\Phi_1$& $\exp\left\{\beta_2^{\top} x_2 + \frac{\sigma^2}{2}\right\}$\\
100d &$\blacksquare$&$\square$&$\square$&$\blacksquare$&& & $1 - \Phi_1$& $\frac{1}{\sigma y}\phi\left(\frac{\ln y -\beta_2^{\top}x_2}{\sigma}\right)\Phi\left(\frac{\beta_1^{\top}x_1+\rho_{12}\frac{\ln y - \beta_2^{\top}x_2}{\sigma}}{\sqrt{1-\rho_{12}^2}}\right)$ & $\exp\left\{\beta_2^{\top} x_2 + \frac{\sigma^2}{2}\right\}\frac{\Phi\left(\beta_1^{\top}x_1+\sigma\rho_{12}\right)}{\Phi_1}$\\
010 &$\square$&$\blacksquare$&$\square$&&&& $1 - \Phi_2$ & $\frac{1}{\sigma}\phi\left(\frac{y-\beta_2^{\top}x_2}{\sigma}\right)$ & $\beta_2^{\top}x_2+\sigma\frac{\phi_2}{\Phi_2}$\\
001i &$\square$&$\square$&$\blacksquare$&&&$\square$& $1 - \Phi_3$ & $\frac{1}{\sigma y}\phi\left(\frac{\ln(\Phi_3 y)-\beta_2^{\top}x_2}{\sigma}\right)\Phi_3$ & $\exp\left\{\beta_2^{\top} x_2 + \frac{\sigma^2}{2}\right\}\frac{1}{\Phi_3}$ \\
001d &$\square$&$\square$&$\blacksquare$&&&$\blacksquare$& $1 - \Phi_3$ & $\frac{1}{\sigma y}\phi\left(\frac{\ln (\Phi_3 y)-\beta_2^{\top}x_2}{\sigma}\right)\Phi\left(\frac{\beta_3^{\top}x_3 + \rho_{23}\frac{\ln (\Phi_3 y) - \beta_2^{\top}x_2}{\sigma}}{\sqrt{1-\rho_{23}^2}}\right)$  & $\exp\left\{\beta_2^{\top} x_2
+ \frac{\sigma^2}{2}\right\}\frac{\Phi\left(\beta_3^{\top}x_3+\sigma\rho_{23}\right)}{\Phi_3^2}$\\
110i &$\blacksquare$&$\blacksquare$&$\square$&$\square$&&& $1 - \Phi_1\Phi_2$ & $\frac{1}{\sigma}\phi\left(\frac{y-\beta_2^{\top}x_2}{\sigma}\right)\Phi_1$ & $\beta_2^{\top} x_2 + \sigma\frac{\phi_2}{\Phi_{2}}$\\
110d &$\blacksquare$&$\blacksquare$&$\square$&$\blacksquare$&&&$1 - \Phi_{12}$ &
$\frac{1}{\sigma}\phi\left(\frac{y-\beta_2^{\top}x_2}{\sigma}\right)\Phi\left(\frac{\beta_1^{\top}x_1 + \rho_{12}\frac{y-\beta_2^{\top}x_2}{\sigma}}{\sqrt{1-\rho_{12}^2}}\right)$ & $\beta_2^{\top} x_2 + \sigma \frac{\Psi_{2|1}}{\Phi_{12}}$\\
101iii &$\blacksquare$&$\square$&$\blacksquare$&$\square$&$\square$&$\square$& $1 - \Phi_1 \Phi_3$ & $\frac{1}{\sigma y}\phi\left(\frac{\ln (\Phi_3 y) -\beta_2^{\top}x_2}{\sigma}\right)\Phi_1\Phi_3$ &$\exp\left\{\beta_2^{\top}x_2+\frac{\sigma^2}{2}\right\}\frac{1}{\Phi_3}$ \\
101dii &$\blacksquare$&$\square$&$\blacksquare$&$\blacksquare$&$\square$&$\square$& $1 - \Phi_1 \Phi_3$ &  $\frac{1}{\sigma y}\phi\left(\frac{(\ln \Phi_3 y) -\beta_2^{\top}x_2}{\sigma}\right)\Phi\left(\frac{\beta_1^{\top}x_1+\rho_{12}\frac{(\ln y \Phi_3) - \beta_2^{\top}x_2}{\sigma}}{\sqrt{1-\rho_{12}^2}}\right)\Phi_3$&
$\exp\left\{\beta_2^{\top} x_2 + \frac{\sigma^2}{2}\right\}\frac{\Phi\left(\beta_1^{\top}x_1+\sigma\rho_{12}\right)}{\Phi_1\Phi_3}$\\
101iid &$\blacksquare$&$\square$&$\blacksquare$&$\square$&$\square$ &$\blacksquare$& $1 - \Phi_1 \Phi_3$ & $\frac{1}{\sigma y}\phi\left(\frac{\ln (\Phi_3 y)-\beta_2^{\top}x_2}{\sigma}\right)\Phi_1\Phi\left(\frac{\beta_3^{\top}x_3 + \rho_{23}\frac{\ln (\Phi_3 y) - \beta_2^{\top}x_2}{\sigma}}{\sqrt{1-\rho_{23}^2}}\right)$& $\exp\left\{\beta_2^{\top} x_2
+\frac{\sigma^2}{2}\right\}\frac{\Phi\left(\beta_3^{\top}x_3+\sigma\rho_{23}\right)}{\Phi_3^2}$\\
011i &$\square$&$\blacksquare$&$\blacksquare$&&&$\square$& $1 - \Phi_2 \Phi_3$& $\frac{1}{\sigma}\phi\left(\frac{\Phi_3y-\beta_2^{\top}x_2}{\sigma}\right)\Phi_3^2$ & $\frac{\beta_2^{\top} x_2}{\Phi_3}
+\sigma \frac{\phi_2}{\Phi_2 \Phi_3}$\\
011d &$\square$&$\blacksquare$&$\blacksquare$&&&$\blacksquare$& $1 - \Phi_{23}$& $\frac{1}{\sigma}\phi\left(\frac{\Phi_3 y-\beta_2^{\top}x_2}{\sigma}\right)\Phi\left(\frac{\beta_3^{\top}x_3 + \rho_{23}\frac{\Phi_3 y - \beta_2^{\top}x_2}{\sigma}}{\sqrt{1-\rho_{23}^2}}\right)\Phi_3$ & $\frac{\beta_2^{\top} x_2}{\Phi_3}
+ \sigma \frac{\Psi_{2|3}}{\Phi_{23}\Phi_3}$ \\
111iii  &$\blacksquare$&$\blacksquare$&$\blacksquare$&$\square$&$\square$&$\square$& $1-\Phi_1\Phi_2\Phi_3$& $\frac{1}{\sigma}\phi\left(\frac{\Phi_3 y-\beta_2^{\top}x_2}{\sigma}\right)\Phi_1\Phi_3^2$& $\frac{\beta_2^{\top} x_2}{\Phi_3} + \sigma \frac{\phi_2}{\Phi_2 \Phi_3}$ \\
111dii &$\blacksquare$&$\blacksquare$&$\blacksquare$&$\blacksquare$&$\square$&$\square$&$1 - \Phi_{12}\Phi_3$& $\frac{1}{\sigma}\phi\left(\frac{\Phi_3 y-\beta_2^{\top}x_2}{\sigma}\right)\Phi\left(\frac{\beta_1^{\top}x_1 + \rho_{12}\frac{\Phi_3 y - \beta_2^{\top}x_2}{\sigma}}{\sqrt{1-\rho_{12}^2}}\right)\Phi_3^2$ & $\frac{\beta_2^{\top} x_2}{\Phi_3}
+ \sigma \frac{\Psi_{2|1}}{\Phi_{12} \Phi_3}$ \\
111iid &$\blacksquare$&$\blacksquare$&$\blacksquare$&$\square$&$\square$ &$\blacksquare$&$1 - \Phi_1\Phi_{23}$& $\frac{1}{\sigma}\phi\left(\frac{\Phi_3 y-\beta_2^{\top}x_2}{\sigma}\right)\Phi_1\Phi\left(\frac{\beta_3^{\top}x_3 + \rho_{23}\frac{\Phi_3 y - \beta_2^{\top}x_2}{\sigma}}{\sqrt{1-\rho_{23}^2}}\right)\Phi_3$ & $\frac{\beta_2^{\top} x_2}{\Phi_3}
+ \sigma \frac{\Psi_{2|3}}{\Phi_{23} \Phi_3}$  \\ \hline\hline
\end{tabular}

\vspace{.5cm}

A blackened square indicates which hurdle or correlation is assumed to
be at work in the model ; an empty square indicates a hurdle which is
not in effect or a zero correlation.

\end{sidewaystable}

The density function $f_+(y)$ may be computed by performing: first the
change of variable $y_2^* = P(I_3=1)y = \Phi_3y$ on the joint density
function of the latent variables entering the hurdle model; then by
integrating this transformed density function over the positive values
of latent variables $y_1^*$ and $y_3^*$.

\begin{itemize}
\item For trivariate hurdle model \code{111} this transformed density function
  is written as:

  \begin{equation}\label{eq13}
    \frac{\Phi_3}{\sigma}\phi\left(y_1^*-\beta_1^\top x_1,\frac{\Phi_3 y-\beta_2^\top x_2}
      {\sigma},y_3^*-\beta_3^\top x_3;\rho_{12},\rho_{13},\rho_{23}\right).
  \end{equation}

  To perform the analytical integration of this function, it is useful
  to rewrite it as the product of the marginal distribution of $y$,
  namely:

  
  \begin{equation}\label{eq14}
    \frac{\Phi_3}{\sigma}\phi\left(\frac{\Phi_3 y-\beta_2^\top x_2}{\sigma}\right)
  \end{equation}

  and of the joint density function of $y_1^*$ and $y_3^*$ conditioned
  with respect to $y$, which can be written as follows:

  \begin{equation}\label{eq15}
    \frac{1}{\sigma_{1|2}\sigma_{3|2}}\phi\left(\frac{y_1^*-\mu_{1|2}}{\sigma_{1|2}},
      \frac{y_3^*-\mu_{3|2}}{\sigma_{3|2}}; \rho_{13|2}\right),
  \end{equation}

  with:
  $$
  \mu_{1|2}=\beta_1^\top x_1+\rho_{12}\frac{\Phi_3y-\beta_2^\top
    x_2}{\sigma}, \quad \mu_{3|2}=\beta_3^\top
  x_3+\rho_{23}\frac{\Phi_3y-\beta_2^\top x_2}{\sigma},
  $$
  $$
  \sigma_{1|2}^2=1-\rho_{12}^2, \quad \sigma_{3|2}^2=1-\rho_{23}^2,
  \quad
  \rho_{13|2}=\frac{\rho_{13}-\rho_{12}\rho_{23}}{\sqrt{1-\rho_{12}^2}\sqrt{1-\rho_{23}^2}}.
  $$

  Using this factorization of the density function of $y_1^*$, $y$ and
  $y_3^*$, we obtain:


  \begin{equation}\label{eq16}
  \begin{array}{l}
    f_+(y)={\displaystyle\frac{\Phi_3}{\sigma}\phi\left(\frac{\Phi_3 y-\beta_2^\top x_2}{\sigma}\right)}\\[12pt]
    \phantom{f_+(y)} \times {\displaystyle\int_0^\infty\int_0^\infty \frac{1}{\sigma_{1|2}\sigma_{3|2}}\phi\left(\frac{y_1^*-\mu_{1|2}}{\sigma_{1|2}},
        \frac{y_3^*-\mu_{3|2}}{\sigma_{3|2}}; \rho_{13|2}\right)dy_1^*dy_3^*} \\[16pt]
    \phantom{f_+(y)}= {\displaystyle\frac{\Phi_3}{\sigma}\phi\left(\frac{\Phi_3 y-\beta_2^\top x_2}{\sigma}\right)
      \int_{-\frac{\mu_{1|2}}{\sigma_{1|2}}}^\infty\int_{-\frac{\mu_{3|2}}{\sigma_{3|2}}}^\infty
      \phi(z_1,z_3;\rho_{13|2})dz_1 dz_3} \\[8pt]
    \phantom{f_+(y)}={\displaystyle\frac{\Phi_3}{\sigma}\phi\left(\frac{\Phi_3 y-\beta_2^\top x_2}{\sigma}\right)}\\[12pt]
    \phantom{f_+(y)} \times{\displaystyle \Phi\left(\frac{\beta_1^\top x_1+\rho_{12}\frac{\Phi_3y-\beta_2^\top x_2}{\sigma}}
        {\sqrt{1-\rho_{12}^2}},\frac{\beta_3^\top x_3+\rho_{23}\frac{\Phi_3y-\beta_2^\top x_2}{\sigma}}
        {\sqrt{1-\rho_{23}^2}};\rho_{13|2}\right)}.
  \end{array}
  \end{equation}


\item For trivariate hurdle model \code{101}, we proceed as for hurdle
  model \code{111} by substituting the joint normal density function
  (\ref{eq13}), by the following joint normal/log-normal density
  function:

  \begin{equation}\label{eq17}
  \frac{1}{\sigma y}\phi\left(y_1^*-\beta_1^\top x_1,\frac{\ln(\Phi_3 y)-\beta_2^\top x_2}
  {\sigma},y_3^*-\beta_3^\top x_3;\rho_{12},\rho_{13},\rho_{23}\right).
  \end{equation}

  To integrate this density function with respect to the positive values of
  $y_1^*$ and $y_2^*$, we rewrite it as the product of the marginal distribution
  of $y$, which is log-normal:

  \begin{equation}\label{eq18}
  \frac{1}{\sigma y}\phi\left(\frac{\ln(\Phi_3 y)-\beta_2^\top x_2}{\sigma}\right)
  \end{equation}

  and of the joint density function of $y_1^*|y$ and $y_3^*|y$, which is bivariate normal:

  \begin{equation}\label{eq19}
  \frac{1}{\sigma_{1|2}\sigma_{3|2}}\phi\left(\frac{y_1^*-\mu_{1|2}}{\sigma_{1|2}},
  \frac{y_3^*-\mu_{3|2}}{\sigma_{3|2}}; \rho_{13|2}\right),
  \end{equation}
  with:
  $$
  \mu_{1|2}=\beta_1^\top x_1+\rho_{12}\frac{\ln(\Phi_3 y)-\beta_2^\top x_2}{\sigma}, \quad
  \mu_{3|2}=\beta_3^\top x_3+\rho_{23}\frac{\ln(\Phi_3 y)-\beta_2^\top x_2}{\sigma},
  $$
  $$
  \sigma_{1|2}^2=1-\rho_{12}^2, \quad \sigma_{3|2}^2=1-\rho_{23}^2, \quad
  \rho_{13|2}=\frac{\rho_{13}-\rho_{12}\rho_{23}}{\sqrt{1-\rho_{12}^2}\sqrt{1-\rho_{23}^2}}.
  $$

  By integrating this factorisation of the density function of $y_1^*$, $y$ and $y_3^*$,
  over the positive values of $y_1^*$ and $y_3^*$, we obtain:

  \begin{equation}\label{eq20}
  \begin{array}{l}
    f_+(y)={\displaystyle\frac{\phi\left(\frac{\ln(\Phi_3 y)-\beta_2^\top x_2}{\sigma}\right)}{\sigma y}
      \int_{-\frac{\mu_{1|2}}{\sigma_{1|2}}}^\infty\int_{-\frac{\mu_{3|2}}{\sigma_{3|2}}}^\infty
      \phi(z_1,z_3;\rho_{13|2})dz_1 dz_3} \\[8pt]
    \phantom{f_+(y)}={\displaystyle\frac{\phi\left(\frac{\ln(\Phi_3 y)-\beta_2^\top x_2}{\sigma}\right)}{\sigma y}}\\[8pt]
    \phantom{f_+(y)}\times{\displaystyle\Phi\left(\frac{\beta_1^\top x_1+\rho_{12}\frac{\ln(\Phi_3 y)-\beta_2^\top x_2}{\sigma}}
        {\sqrt{1-\rho_{12}^2}},\frac{\beta_3^\top x_3+\rho_{23}\frac{\ln(\Phi_3 y)-\beta_2^\top x_2}{\sigma}}
        {\sqrt{1-\rho_{23}^2}};\rho_{13|2}\right)}.
  \end{array}
  \end{equation}

\item The density function $f_+(y)$ for bivariate hurdle models
  \code{011} and \code{110} and univariate hurdle model \code{010} can
  be derived from that of trivariate model \code{111} by eliminating
  hurdles 1, 3, 1 and 3, respectively.  Likewise, this density
  function for bivariate hurdle models \code{100} and \code{001} can
  be derived from that of trivariate hurdle model \code{101} by
  eliminating hurdles 1 and 3, respectively.  Corresponding formulas
  for $f_+(y)$ for all this special cases implemented in \proglang{R}
  are presented in Table \ref{tab:typo}.
\end{itemize}

From these results it is easy to derive the likelihood function of a
random sample of $n$ observations of the censored dependent variable
$y$.  As these observations are all independently drawn from the same
conditional (on covariates $x_1$, $x_2$ and $x_3$) discrete-continuous
distribution, which assigns a conditional probability mass $P(y=0)$ to
the observed value $y=0$ and a conditional density function $f_+(y)$
to the observed values $y>0$, the log-likelihood function for an
observation $y_i$ can be written as :

\begin{equation}\label{eq21}
\ln L_i = \left\{
  \begin{array}{ll}
    \ln P(y_i=0) & \mbox{if} \quad y_i=0 \\
    \ln f_+(y_i) & \mbox{if} \quad y_i>0
\end{array}
\right.
\end{equation}

and the log-likelihood for the entire random sample:

\begin{equation}\label{eq22}
  \ln L = \sum_{i=1}^n \ln L_i= \sum_{i \mid y_i = 0} \ln P(y_i=0) + \sum_{i \mid y_i > 0} \ln f_+(y_i).
\end{equation}

\subsection{Heteroscedasticity and nonnormality}

\change{
Contrary to the classical linear regression model which estimation 
is robust with respect to deviations from the assumptions of homoscedasticity
and normality of disturbances, the maximum likelihood estimation of a Tobit model
become inconsistent under heteroscedasticity and nonnormality of disturbances.
Therefore, it is important to have methods allowing to test whether these 
assumptions are acceptable, on grounds of empirical evidence provided by 
a sample of observations, and to suggest how to respecify the model in case where
a misspecification is brought out. In this section, we shall tackle this problem by
using more flexible specifications of the desired consumption relation, where 
homoscedasticity and normality assumptions appear as special cases. Our choices are
inspired by model specifications identified in our survey of hurdle Tobit model 
applications reported in Table~\ref{tab:exconso}.
}
Stated in general terms, our generalizations of models presented in section 2.1
are written as:

\begin{equation}\label{eq51}
  \left\{
  \begin{array}{l}
  y_1^* = \beta_1^\top x_1 + z_1 \\[4pt]
  T(y_2^*) = \beta_2^\top x_2 + \sigma(\beta_0^\top x_0) z_2 \\[4pt]
  y_3^* = \beta_3^\top x_3 + z_3 \\
  \end{array}
  \right.
\end{equation}

where $T(y_2^*)$ denotes a monotonic transformation of desired consumption 
variable $y_2^*$, $\sigma(\beta_0^\top x_0)$ a positive monotonic transformation 
of a linear function of a vector of covariates $x_0$, selected to explain
the heteroscedasticity of desired consumption, $z_1$, $z_2$ and $z_3$ 
standard normal random variables, with $z_2$ possibly truncated at the bounds of an interval
$]B_1,B_2[$ ensuring that the range of values of $\beta_2^\top x_2 + \sigma(\beta_0^\top x_0) z_2$
corresponds to the domain of definition of the inverse transformation
$y_2^* = T^{-1} (\beta_2^\top x_2 + \sigma(\beta_0^\top x_0) z_2)$.

With these assumptions in mind, we can derive the form of the observable dependent variable

\begin{equation}\label{eq52}
  y = \frac{T^{-1} (\beta_2^\top x_2 + \sigma(\beta_0^\top x_0) z_2)}{P(I_3=1)} I_1 I_2 I_3 .
\end{equation}

from the joint distribution of random variables $z_1$, $z_2$ and $z_3$, which is written as:

\begin{equation}\label{eq53}
    \frac{\phi\left(z_1,z_2,z_3;\rho_{12},\rho_{13},\rho_{23}\right)}{\Phi(B_2)-\Phi(B_1)}.
\end{equation}

Using the change of variables:

\begin{equation}\label{eq54}
  \left\{
  \begin{array}{l}
   z_1  =  y_1^* - \beta_1^\top x_1 \\[4pt]
   z_2  =  {\displaystyle\frac{T(y_2^*) - \beta_2^\top x_2}{\sigma(\beta_0^\top x_0)}} \\[8pt]
   z_3  =  y_3^* - \beta_3^\top x_3 \\
  \end{array}
  \right.
\end{equation}

we derive the joint density function of latent variables $y_1^*$, $y_2^*$ and $y_3^*$:

\begin{equation}\label{eq55}
  f(y_1^*, y_2^*, y_3^*)=\frac{T^{'}(y_2^*)}{\sigma(\beta_0^\top x_0)}
  \frac{\phi\left(y_1^* - \beta_1^\top x_1, \frac{T(y_2^*)-\beta_2^\top x_2}{\sigma(\beta_0^\top x_0)}, 
  y_3^* - \beta_3^\top x_3)\right)}{\Phi(B_2)-\Phi(B_1)},
\end{equation}

where $T^{'}(y_2^*)$ stands for the derivative of $T(y_2^*)$. Using this density function 
we compute the probability mass $P(y=0)=1-P(y>0)$ as follows:

\begin{equation}\label{eq56}
  \begin{array}{l}
  P(y=0) = {\displaystyle 1-\int_{0}^\infty \int_{0}^\infty \int_{0}^\infty \frac{\phi(y_1^*,y_2^*,y_3^*;\rho_{12},\rho_{13},\rho_{23})}
  {\Phi(B_2)-\Phi(B_1)} dy_1^* dy_2^* dy_3^*}\\[16pt]
  \phantom{P(y=0)} ={\displaystyle 1-\int_{-\beta_1^\top x_1}^\infty \int_{\frac{T(0)-\beta_2^\top x_2}{\sigma(\beta_0^\top x_0)}}^{B_2}
  \int_{-\beta_3^\top x_3}^\infty \frac{\phi(z_1,z_2,z_3;\rho_{12},\rho_{13},\rho_{23})}{\Phi(B_2)-\Phi(B_1)} dz_1 dz_2 dz_3}\\[16pt]
  \phantom{P(y=0)} ={\displaystyle 1-\frac{\Phi(\beta_1^\top x_1,\frac{\beta_2^\top x_2 - T(0)}{\sigma(\beta_0^\top x_0)},
  \beta_3^\top x_3;\rho_{12},\rho_{13},\rho_{23})- \Phi(\beta_1^\top x_1,-B_2,\beta_3^\top x_3;\rho_{12},\rho_{13},\rho_{23})}
  {\Phi(B_2)-\Phi(B_1)}}.
  \end{array}
  \end{equation}

To compute the density function $f_+(y)$ , we first perform the
change of variable $y_2^* = \Phi_3y$ on the joint density
function \label{eq55}, leading to the joint density function of variables $y_1^*$, $y$ and $y_3^*$:

\begin{equation}\label{eq57}
  f(y_1^*, y, y_3^*)=\frac{\Phi_3 T^{'}(\Phi_3 y)}{\sigma(\beta_0^\top x_0)}
  \frac{\phi\left(y_1^* - \beta_1^\top x_1, \frac{T(\Phi_3 y)-\beta_2^\top x_2}{\sigma(\beta_0^\top x_0)}, 
  y_3^* - \beta_3^\top x_3)\right)}{\Phi(B_2)-\Phi(B_1)}.
\end{equation}

Then we integrate this transformed density function, rewritten as the product 
of the marginal density function of $y$, namely:

\begin{equation}\label{eq58}
    \frac{\Phi_3 T^{'}(\Phi_3 y)}{\sigma(\beta_0^\top x_0)} 
    \frac{\phi\left(\frac{T(\Phi_3 y)-\beta_2^\top x_2}{\sigma(\beta_0^\top x_0)}\right)}{\Phi(B_2)-\Phi(B_1)},
\end{equation}

and of the joint density function of $y_1^*$ and $y_2^*$ conditioned with respect to $y$,
stated by formula \label{eq15}, over the positive values of latent variables $y_1^*$ and $y_3^*$.
By this way we get:


\begin{equation}\label{eq59}
  \begin{array}{l}
    f_+(y)={\displaystyle 
    \frac{\Phi_3 T^{'}(\Phi_3 y)}{\sigma(\beta_0^\top x_0)} 
    \frac{\phi\left(\frac{T(\Phi_3 y)-\beta_2^\top x_2}{\sigma(\beta_0^\top x_0)}\right)}{\Phi(B_2)-\Phi(B_1)}}\\%[12pt]
    \phantom{f_+(y)} \times{\displaystyle \Phi\left(\frac{\beta_1^\top x_1+\rho_{12}\frac{T(\Phi_3 y)-\beta_2^\top x_2}
    {\sigma(\beta_0^\top x_0)}}{\sqrt{1-\rho_{12}^2}},\frac{\beta_3^\top x_3+\rho_{23}\frac{T(\Phi_3 y)-\beta_2^\top x_2}
    {\sigma(\beta_0^\top x_0)}}{\sqrt{1-\rho_{23}^2}};\rho_{13|2}\right)}.
  \end{array}
\end{equation}

A natural choice for the heteroscedastic model $\sigma(\beta_0^\top x_0)$ is given by the exponential 
functional form $\sigma(\beta_0^\top x_0)=\exp\{\beta_0^\top x_0\}$, as exponential is a 
strictly increasing function mapping the set of real numbers $]-\infty,+\infty[$ onto the
the set of positive real numbers $]0,+\infty[$, and leads to homoscedasticity when the linear function
$\beta_0^\top x_0$ is a constant. This allows to test the assumption of homoscedasticity by inserting
an intercept among parameters $\beta_0$ and assessing whether the other parameters are statistically 
significant or not. 

The main drawbacks of this functional form lies in its implicit assumption of
unboundedness of the variance of disturbance $\epsilon_2$ with respect to covariates $x_0$. When a priori
information suggests that this conditional variance remains bounded with respect to any possible values of
these covariates, it is best to consider a functional model of the form 
$\sigma(\beta_0^\top x_0)=\exp\{\alpha F(\beta_0^\top x_0)\}$, where $F(.)$ is a distribution function mapping 
the set of real numbers $]-\infty,+\infty[$ onto the unit interval $]0,1[$, and $\alpha$ a parameter playing 
the role of the intercept in the exponential functional form. Therefore, testing the assumption of homoscedasticity 
amounts to assessing the statistical non significance of parameter vector $\beta_0$ without an intercept, and setting 
the standard deviation of $\epsilon_2$ equal to $\sigma = \exp\{\alpha F(0)\}$. For a practical application of this
parametric model of heteroscedasticity, a menu of functional forms of $F(.)$ have been programmed in \pkg{mhurdle},
including the distribution functions of standard normal, logistic, Cauchy and Gompertz random variables.

As far as the choice of transformation $T(y_2^*)$ is concerned, two families of parametric functions have been
considered, in order to generate departures from normality towards skewed and leptokurtic (more sharply peaked)
distributions, of the kind encountered in collected economic data.

To generate skewed distributions of $y_2^*$ we use the two parameters \citet{BOX:COX:64} transformation, as suggested
by \citet{CHAZ:05}. This transformation is written as:

\begin{equation}\label{eq60}
  T(y_2^*) = \left\{
  \begin{array}{ll}
    \frac{(y_2^* + \gamma)^{\lambda}-1}{\lambda} & \mbox{if} \quad \lambda \neq 0 \\
    \ln (y_2^* + \gamma) & \mbox{if} \quad \lambda=0
  \end{array}
  \right.
\end{equation}

with $\lambda$ a parameter characterizing the non linearity of the transformation. 
As shown in Figure 1, this transformation is linear for $\lambda = 1$, convex for $\lambda > 1$, and concave for
$\lambda < 1$, with a ceiling asymptote at $-1/\lambda$ when $\lambda<0$. $\gamma$ is a 
location parameter, restricting the domain where the transformation holds, to the 
interval $]-\lambda,+\infty[$. Hence, the image of this interval by the Box-Cox 
transformation is given by:

\begin{equation}\label{eq61}
  T(]-\gamma,+\infty[) = \left\{
  \begin{array}{ll}
    ]-1/\lambda,+\infty[ & \mbox{if} \quad \lambda > 0 \\
    ]-\infty,+\infty[ & \mbox{if} \quad \lambda = 0 \\
    ]-\infty,-1/\lambda[ & \mbox{if} \quad \lambda <0 
  \end{array}
  \right.
\end{equation}

As the inverse Box-Cox transformation is written as:

\begin{equation}\label{eq62}
  y_2^* = \left\{
  \begin{array}{ll}
    \{\lambda(\beta_2^\top x_2 + \sigma(\beta_0^\top x_0) z_2)+1\}^{1/\lambda}-\gamma & \mbox{if} \quad \lambda \neq 0 \\
    \exp\{\beta_2^\top x_2 + \sigma(\beta_0^\top x_0) z_2\} - \gamma & \mbox{if} \quad \lambda=0
  \end{array}
  \right.
\end{equation}

it turns out that $z_2$ is truncated at the bounds of the interval $]B_1,B_2[$, with:

\begin{equation}\label{eq63}
  B_1 = \left\{
  \begin{array}{ll}
    -\frac{\beta_1^\top x_1 + (1/\lambda)}{\sigma(\beta_0^\top x_0)} & \mbox{if} \quad \lambda > 0 \\
    -\infty & \mbox{if} \quad \lambda \leq 0
  \end{array}
  \right.
\end{equation}

and

\begin{equation}\label{eq64}
  B_2 = \left\{
  \begin{array}{ll}
     +\infty & \mbox{if} \quad \lambda \geq 0\\
    -\frac{\beta_1^\top x_1 + (1/\lambda)}{\sigma(\beta_0^\top x_0)} & \mbox{if} \quad \lambda < 0
  \end{array}
  \right.
\end{equation}

Hence:

\begin{equation}\label{eq63}
  \Phi(B_2)-\Phi(B_1) = \left\{
  \begin{array}{ll}
    \Phi(sign(\lambda)\frac{\beta_1^\top x_1 + (1/\lambda)}{\sigma(\beta_0^\top x_0)}) & \mbox{if} \quad \lambda \neq 0 \\
    1 & \mbox{if} \quad \lambda = 0
  \end{array}
  \right.
\end{equation}

with

\begin{equation}\label{eq69}
  sign(\lambda) = \left\{
  \begin{array}{ll}
     + & \mbox{if} \quad \lambda > 0\\
     - & \mbox{if} \quad \lambda < 0
  \end{array}
  \right.
\end{equation}
 
Finally, the distribution of the observed variable $y$ is obtained by inserting these results 
in formulas \label{eq56} and \label{eq59}, jointly with the derivative of the Box-Cox 
transformation, $T^{'}(y_2^*)=(y_2^*+\gamma)^{\lambda-1}$.

Note that when $\gamma = 0$, the Box-Cox transformation holds only for $y_2^* > 0$, corresponding 
to the restriction enforced to the desired consumption relation, when hurdle 2 is not in effect, 
by means of Cragg's log-normal or truncated normal functional forms. These functional forms arise as special cases
of the Box-Cox transformation, when $\lambda = 0$ and $\lambda = 1$, respectively. Thus, testing the
statistical significance of parameter $\gamma$, against the alternative $\gamma > 0$, amounts
to testing the assumption that hurdle 2 is not in effect within the framework of a nested
model. Conversely, when $\gamma < 0$ the Box-Cox transformation holds for $y_2^* > -\gamma$
meaning that $-\gamma$ stands for a "committed" consumption of a basic necessity, while 
$\gamma > 0$ typify a luxury good whose consumption occurs only above a given income threshold.

The way the Box-Cox transformation generates the profile of the density function of $y_2^*$ from that
of $z_2$ can be illustrated by considering the special case of $\beta_2^\top x_2=0$ and
$\sigma(\beta_0^\top x_0)=1$. In this special case, the marginal density function of $y_2^*$
is defined for $y_2^* \geq 0$,and takes the form of the product of the density function of $z_2$,
written as a function of $y_2^*$:

\begin{equation}\label{eq64}
   \frac{\phi(\frac{y_2^{*\lambda}-1}{\lambda})}{\Phi(sign(\lambda)\frac{1}{\lambda})}
\end{equation}

times the value of the Jacobian of the Box-Cox transformation:

\begin{equation}\label{eq65}
   J = y_2^{*(\lambda - 1)}.
\end{equation}

This Jacobian magnifies or reduces the value of the density function of $z_2$ depending on 
whether $J>1$ or $J<1$, respectively. According to the values of $\lambda$, three cases 
must be considered.

\begin{itemize}

\item When $\lambda = 1$, $J=1$ whatever the value of $y_2^*$. Hence, the shape of the density 
 function of $y_2^*$ is the same as that of $z_2$, namely a normal random variable left truncated
 at $y_2^* = 0$ and with a mode at $y_2^* = 1$.

\item When $\lambda > 1$, $J$ is an unbounded increasing function of $y_2^*$, taking value $0$ at 
 $y_2^* = 0$ and $1$ at $y_2^* = 1$, while the density function of $z_2$ at $y_2^*=0$ holds finite. 
 Hence, the density function of $y_2^*$ starts from $0$ and shifts the mode beyond $y_2^* = 1$. 
 Furthermore, the Jacobian being a linear function of
 $y_2^*$ when $\lambda = 2$, a concave function when $1 < \lambda < 2$, and a convex function 
 when $\lambda > 2$, the concentration of the density function around its mode increases with
 the value of $\lambda$, until the collapse of the entire probability mass at $y_2^* = 1$, when
 $\lambda \rightarrow \infty$.    

\item When $\lambda < 1$, $J$ is a decreasing function of $y_2^*$, from $+\infty$ at $y_2^* = 0$,
 to $0$ when $y_2^* \rightarrow \infty$, taking value $1$ at $y_2^* = 1$, while the density function 
 of $z_2$ at $y_2^*=0$ holds finite or takes value $0$ depending on whether $0<\lambda<1$ or $\lambda\leq 0$. 
 Hence, contrary to the previous case, the mode of the density function of $y_2^*$ shifts short of $y_2^* = 1$, 
 while, at the origin, the density function of $y_2^*$ tends to $+\infty$ and turns out to be an 
 undetermined product depending on whether $0<\lambda<1$ or $\lambda\leq 0$. The application of
 Hospital's rule allows to remove the indetermination of this limit, which turns out to be equal to $0$.
 Note that the same tendency towards a concentration of the density function of $y_2^*$ around its degenerated 
 limit at $y_2^* = 1$ appears when $\lambda \rightarrow -\infty$. Note also that the J-shaped profile
 of this density function, when $0<\lambda<1$, has a mode preceded by an antimode that can coalesce into 
 a point of inflexion.
 
\end{itemize}
 
All these profiles are illustrated in Figures 2.

To generate leptokurtic distributions of $y_2^*$ we use the transformation popularized by \citet{JOHN:49},
based on the inverse hyperbolic sine function, namely:

\begin{equation}\label{eq66}
  T(y_2^*) = \frac{1}{\theta} \sinh^{-1}(\theta y_2^*) = \frac{1}{\theta} 
  \ln \{\theta y_2^* + \sqrt{(\theta y_2^*)^2+1}\} 
\end{equation}

with $\theta$ a parameter characterizing the non linearity of the transformation. Indeed, as illustrated
by Figure 3, while, for $\theta > 0$ and finite, the general shape of this transformation is that of 
an odd \footnote{A function $f(x)$ is said odd if $f(-x)=-f(x)$ whatever $x$, implying that $f(0)=0$.} increasing 
function, from $-\infty$ to $+\infty$,  with an inflexion point at $y_2^*=0$, concave for $y_2^*>0$ and 
convex for $y_2^*<0$, the transformation becomes linear, 
when $\theta \rightarrow 0$, and constant (equal to $0$), when $\theta\rightarrow +\infty$ 
\footnote{These limits can be obtained easily by using Hospital's rule.}, by progressively loosing its 
sigmoidal profile. As these profiles are invariant with respect to a change of sign of $\theta$, this parameter
can be assumed to be non negative, without loss of generality.

To determine the distribution of the observed variable $y$, we need to insert in formulas \label{eq56} and \label{eq59}
the expressions of $\Phi(B_2)-\Phi(B_1)$ and $T^{'}(y_2^*)$ for Johnson's transformation. We first note that 
the inverse function of this transformation, which is written:

\begin{equation}\label{eq67}
  \begin{array}{l}
  y_2^* = {\displaystyle \frac{\sinh (\theta(\beta_2^\top x_2 + \sigma(\beta_0^\top x_0) z_2))}{\theta} }\\[16pt]
  \phantom{ y_2^*} = {\displaystyle \frac{\exp\{\theta(\beta_2^\top x_2 + \sigma(\beta_0^\top x_0) z_2)\}
  -\exp\{-\theta(\beta_2^\top x_2 + \sigma(\beta_0^\top x_0) z_2)\}}{2 \theta}}
  \end{array} 
\end{equation}

does not require a truncation of the support of $z_2$, implying that $B_1=-\infty$, $B_2=+\infty$, and hence
$\Phi(B_2)-\Phi(B_1)=1$. 

Secondly, the derivative of Johnson's transformation, namely:

\begin{equation}\label{eq68}
  T^{'}(y_2^*) = \frac{1}{\sqrt{(\theta y_2^*)^2 + 1}}   
\end{equation}

is a bell-shaped positive pair function \footnote{A function $f(x)$ is said pair if $f(-x)=f(x)$ whatever $x$.} taking
constant values $1$ and $0$ when $\theta \rightarrow 0$ and $\theta\rightarrow +\infty$, respectively. 

To analyze the profile of the density function of $y_2^*$ generated by this transformation,
we consider, as we did before for the Box-Cox transformation, the special case of $\beta_2^\top x_2=0$ and
$\sigma(\beta_0^\top x_0)=1$. In this special case, the marginal density function of $y_2^*$ 
takes the form of the product of the density function of $\phi (T(y_2^*))$ times the value of the Jacobian
$T^{'}(y_2^*)$. As a function of $\theta$, this formula generates a family of bell-shaped density
functions which kurticity increases with the value of $\theta$, from metakurticity (that of the normal distribution) 
when $\theta=0$, to extreme leptokurticity resulting in a concentration of all the density function at point 
$y_2^*=0$, when $\theta \rightarrow +\infty$, as shown in Figure 4. Therefore, the inverse hyperbolic sine 
transformation must be used only when hurdle 2 is in effect.


\section{Model estimation, evaluation and selection}
\label{sec:estevalsel}

The econometric framework described in the previous section provides a
theoretical background for tackling the problems of model estimation,
evaluation and selection within the statistical theory of classical
inference.

\subsection{Model estimation}

The full parametric specification of our multiple hurdle models allows
to efficiently estimate their parameters by means of the maximum
likelihood principle. Indeed, it is well known from classical
estimation theory that, under the assumption of a correct model
specification and for a likelihood function sufficiently well behaved,
the maximum likelihood estimator is asymptotically efficient within
the class of consistent and asymptotically normal estimators
\footnote{See \citet{AMEM/85} chapter 4, for a more rigorous statement
  of this property.}.

More precisely, the asymptotic distribution of the maximum likelihood
estimator $\hat{\theta}$ for the parameter vector $\theta$ of a
multiple hurdle model, is written as:

\begin{equation}\label{eq23}
  \hat{\theta} \overset{A}{\sim} N(\theta,\frac{1}{n} I_A (\theta)^{-1}),
\end{equation}

where $\overset{A}{\sim}$ stands for ``asymptotically distributed as''
and

\begin{equation*}
  I_A(\theta)=\mbox{plim} \frac{1}{n} \sum_{i=1}^n E\left(\frac{\partial^2
      \ln L_i (\theta)}{\partial \theta \partial \theta^\top}\right)
  =\mbox{plim} \frac{1}{n} \sum_{i=1}^n E\left(\frac{\partial\ln
      L_i(\theta)}{\partial\theta} \frac{\partial\ln L_i
      (\theta)}{\partial\theta^\top}\right)
\end{equation*}

for the asymptotic Fisher information matrix of a sample of $n$
independent observations.

More generally, any inference about a differentiable vector function
of $\theta$, denoted by $\gamma=h(\theta)$, can be based on the
asymptotic distribution of its implied maximum likelihood estimator
$\hat{\gamma}=h(\hat{\theta})$.  This distribution can be derived from
the asymptotic distribution of $\hat{\theta}$ according to the so
called delta method:

\begin{equation}\label{eq24}
\hat{\gamma} \overset{A}{\sim} h(\theta)+\frac{\partial
  h}{\partial\theta^\top} (\hat\theta-\theta) \overset{A}{\sim}
N\left(\gamma,\frac{1}{n} \frac{\partial h}{\partial\theta^\top}I_A
(\theta)^{-1}\frac{\partial h^\top}{\partial\theta}\right).
\end{equation}

The practical use of these asymptotic distributions requires to
replace the theoretical variance-covariance matrix of these asymptotic
distributions with consistent estimators, which can be obtained by
using $\frac{\partial h (\hat{\theta})}{\partial\theta^\top}$ as a
consistent estimator for $\frac{\partial h
  (\theta)}{\partial\theta^\top}$ and either $\frac{1}{n} \sum_{i=1}^n
\frac{\partial^2\ln L_i
  (\hat{\theta})}{\partial\theta\partial\theta^\top}$ or $\frac{1}{n}
\sum_{i=1}^n \frac{\partial\ln L_i(\hat{\theta})}{\partial\theta}
\frac{\partial\ln L_i (\hat{\theta})}{\partial\theta^\top}$ as a
consistent estimator for $I_A(\theta)$. The last two estimators are
directly provided by two standard iterative methods used to compute
the maximum likelihood parameter's estimate, namely the Newton-Raphson
method and the Berndt, Hall, Hall, Hausman or \textsc{bhhh} method,
respectively, mentioned in section 4.3.

\subsection{Model evaluation and selection using goodness of fit measures}

Two fundamental principles should be used to appraise the results of a
model estimation, namely its economic relevance and its statistical
and predictive adequacy. The first principle deals with the issues of
accordance of model estimate with the economic rationale underlying
the model specification and of its relevance for answering the
questions for which the model has been built. These issues are
essentially context specific and, therefore, cannot be dealt with by
means of generic criteria.  The second principle refers to the issues
of empirical soundness of model estimate and of its ability to predict
sample or out-of-sample observations.  These issues can be tackled by
means of formal tests of significance, based on the previously
presented asymptotic distributions of model estimates, and by measures
of goodness of fit/prediction, respectively.

To assess the goodness of fit of \pkg{mhurdle} estimates, two pseudo
$R^2$ coefficients are provided. The first one is an extension of the
classical coefficient of determination, used to explain the fraction
of variation of the dependent variable explained by the covariates
included in a linear regression model with intercept. The second one
is an extension of the likelihood ratio index introduced
by \citet{MCFA/74} to measure the relative gain in the maximised
log-likelihood function due to the covariates included in a
qualitative response model.

To define a pseudo coefficient of determination, we rely on the non
linear regression model explaining the dependent variable of a
multiple hurdle model. This model is written as:

\begin{equation}\label{eq25}
y=E(y)+u,
\end{equation}

where $u$ stands for a zero expectation, heteroskedastic random
disturbance and $E(y)$ for the expectation of the censored dependent
variable $y$:

\begin{equation}\label{eq26}
  E(y)=0 \times P(y=0)+\int_0^\infty yf_+(y)dy=\int_0^\infty yf_+(y)dy.
\end{equation}

To compute this expectation, we reformulate it as a multiple integral
of the joint density function of $y_1^*$, $y$ and $y_3^*$ multiplied
by $y$, over the positive values of these variables.

\begin{itemize}

\item For trivariate hurdle model \code{111}, using the density function
   (\ref{eq13}) and the change of variables:

  \begin{equation}\label{eq27}
  \left\{
  \begin{array}{l}
   z_1 = y_1^* - \beta_1^\top x_1 \\[4pt]
   z_2 = {\displaystyle\frac{\Phi_3 y - \beta_2^\top x_2}{\sigma}} \\[8pt]
   z_3 = y_3^* - \beta_3^\top x_3
  \end{array}
  \right.
  \end{equation}

  this reformulation of $E(y)$ is written as:

  \begin{equation}\label{eq28}
  \begin{array}{l}
  E(y)={\displaystyle\int_{-\beta_1^\top x_1}^\infty \int_{-\frac{\beta_2^\top x_2}{\sigma}}^\infty
   \int_{-\beta_3^\top x_3}^\infty \frac{\beta_2^\top x_2+\sigma z_2}{\Phi_3}
   \phi\left(z_1,z_2,z_3;\rho_{12},\rho_{13},\rho_{23}\right) dz_1 dz_2 dz_3} \\[16pt]
 \phantom{E(y)} ={\displaystyle\frac{\beta_2^\top x_2}{\Phi_3}\Phi\left(\beta_1^\top x_1,\frac{\beta_2^\top x_2}{\sigma},
   \beta_3^\top x_3;\rho_{12},\rho_{13},\rho_{23}\right)} \\[12pt]
  \phantom{E(y)}+{\displaystyle\frac{\sigma}{\Phi_3}
   \int_{-\beta_1^\top x_1}^\infty \int_{-\frac{\beta_2^\top x_2}{\sigma}}^\infty
   \int_{-\beta_3^\top x_3}^\infty z_2\phi\left(z_1,z_2,z_3;\rho_{12},\rho_{13},\rho_{23}\right) dz_1 dz_2 dz_3}.
  \end{array}
  \end{equation}

  To perform the analytical integration of the second term of the right-hand side
  of this formula, it is useful to rewrite the density function of $z_1$, $z_2$
  and $z_3$ as the product of the marginal density function of $z_1$ and $z_2$, namely
  $\phi\left(z_1,z_2;\rho_{13}\right)$ and of the density function of $z_2|z_1,z_3$,
  which can be written as follows:

  \begin{equation}\label{eq29}
  \frac{\phi\left(\frac{z_2-\mu_{2|1,3}}{\sigma_{2|1,3}}\right)}{\sigma_{2|1,3}},
  \end{equation}

  where:
  $$
  \mu_{2|1,3}=\varrho_1 z_1 + \varrho_3 z_3, \quad
  \sigma_{2|1,3}^2=\frac{1-\rho_{12}^2-\rho_{13}^2-\rho_{23}^2+2\rho_{12}\rho_{13}\rho_{23}}
  {1-\rho_{13}^2},
  $$
  with:
  $$
  \varrho_1=\frac{\rho_{12}-\rho_{13}\rho_{23}}{1-\rho_{13}^2}, \quad
  \varrho_3=\frac{\rho_{23}-\rho_{12}\rho_{13}}{1-\rho_{13}^2}.
  $$

  Using this factorisation of the density function of $z_1$, $z_2$ and $z_3$,
  we obtain:

  \begin{equation}\label{eq30}
  \begin{array}{l}
  {\displaystyle\int_{-\beta_1^\top x_1}^\infty \int_{-\frac{\beta_2^\top x_2}{\sigma}}^\infty
  \int_{-\beta_3^\top x_3}^\infty z_2\phi\left(z_1,z_2,z_3;\rho_{12},\rho_{13},\rho_{23}\right)
  dz_1 dz_2 dz_3}\\[12pt]
  ={\displaystyle\int_{-\beta_1^\top x_1}^\infty \int_{-\beta_3^\top x_3}^\infty
  \left[\int_{-\frac{\beta_2^\top x_2}{\sigma}}^\infty z_2\phi\left(\frac{z_2-\mu_{2|1,3}}
  {\sigma_{2|1,3}}\right)\frac{dz_2}{\sigma_{2|1,3}}\right]\phi\left(z_1,z_2;\rho_{13}\right)dz_1 dz_3}.
  \end{array}
  \end{equation}

  By performing the change of variable:

  \begin{equation}\label{eq31}
  z=\frac{z_2 - \mu_{2|1,3}}{\sigma_{2|1,3}},
  \end{equation}

  the integral with respect to $z_2$ simplifies to:

  \begin{equation}\label{eq32}
  \mu_{2|1,3}\Phi\left(\frac{\frac{\beta_2^\top x_2}{\sigma}+\mu_{2|1,3}}{\sigma_{2|1,3}}\right)
  +\sigma_{2|1,3}\phi\left(\frac{\frac{\beta_2^\top x_2}{\sigma}+\mu_{2|1,3}}{\sigma_{2|1,3}}\right).
  \end{equation}

  By inserting this result in formula (\ref{eq30}), we finally obtain:

  \begin{equation}\label{eq33}
  \begin{array}{l}
  E(y)=\displaystyle\frac{\beta_2^\top x_2}{\Phi_3 }\Phi\left(\beta_1^\top x_1,\frac{\beta_2^\top x_2}{\sigma},
  \beta_3^\top x_3;\rho_{12},\rho_{13},\rho_{23}\right)\\[8pt]
\displaystyle \phantom{E(y)} + \frac{\sigma}{\Phi_3}
  \int_{-\beta_1^\top x_1}^\infty \int_{-\beta_3^\top x_3}^\infty \left[\left(\varrho_1 z_1+\varrho_3 z_3 \right)
  \Phi\left(\frac{\frac{\beta_2^\top x_2}{\sigma}
  +\varrho_1 z_1+\varrho_3 z_3}{\sigma_{2|1,3}}\right)\right.\\[16pt]
\left.\phantom{E(y)}+\sigma_{2|1,3}
  \phi\left(\frac{\frac{\beta_2^\top x_2}{\sigma}+\varrho_1 z_1+\varrho_3 z_3}{\sigma_{2|1,3}}\right)
  \right]\phi\left(z_1,z_3;\rho_{13}\right) dz_1 dz_3.
  \end{array}
  \end{equation}

\item For trivariate hurdle model \code{101}, we proceed as for hurdle
  model \code{111} by first substituting the joint normal density
  function (\ref{eq13}) by the joint normal/log-normal density
  function (\ref{eq17}), then by performing the change of variables:

  \begin{equation}\label{eq34}
  \left\{
  \begin{array}{l}
   z_1 = y_1^* - \beta_1^\top x_1 \\[4pt]
   z_2 = {\displaystyle\frac{\ln\left(\Phi_3 y\right) - \beta_2^\top x_2}{\sigma}} \\[12pt]
   z_3 = y_3^* - \beta_3^\top x_3
  \end{array}
  \right.
  \end{equation}

  This leads to the following expression of the expected value of $y$:

  \begin{equation}\label{eq35}
  \begin{array}{l}
  E(y)=\displaystyle\int_{-\beta_1^\top x_1}^\infty \int_{-\infty}^\infty
  \int_{-\beta_3^\top x_3}^\infty \frac{\exp\{\beta_2^\top x_2+\sigma z_2\}}{\Phi_3} \\
  \displaystyle  \times   \phi\left(z_1,z_2,z_3;\rho_{12},\rho_{13},\rho_{23}\right) dz_1 dz_2 dz_3
  \displaystyle =\frac{\exp\{\beta_2^\top x_2\}}{\Phi_3}\\[12pt]
  \displaystyle \times \int_{-\beta_1^\top x_1}^\infty \int_{-\beta_3^\top x_3}^\infty 
  \left[\int_{-\infty}^\infty \exp\{\sigma z_2\}\phi\left(\frac{z_2-\mu_{2|1,3}}  
  {\sigma_{2|1,3}}\right)\frac{dz_2}{\sigma_{2|1,3}} \right] \phi\left(z_1,z_3;\rho_{13}\right)dz_1 dz_3
  \end{array}
  \end{equation}

  obtained by factorising the density function of $z_1$, $z_2$ and $z_3$
  as the product of the marginal density function of $z_1$ and $z_3$
  times the density function of $z_2|z_1,z_3$.\\
  By performing the change of variable (\ref{eq31}),
  the integral with respect to $z_2$ simplifies to:

  \begin{equation}\label{eq36}
  \int_{-\infty}^\infty \exp\{\sigma (\mu_{2|1,3}+\sigma_{2|1,3}z)\}\phi(z)dz
  =\exp\left\{\sigma\mu_{2|1,3}+\frac{\sigma^2\sigma_{2|1,3}^2}{2}\right\}.
    \end{equation}

  By inserting this result in formula (\ref{eq35}), we finally obtain:

  \begin{equation}\label{eq37}
 \begin{array}{l}
\displaystyle  E(y)=\frac{\exp\left\{\beta_2^\top x_2+\frac{\sigma^2\sigma_{2|1,3}^2}{2}\right\}}
  {\Phi_3}\\[12pt]
\displaystyle \phantom{E(y)} \times \int_{-\beta_1^\top x_1}^\infty \int_{-\beta_3^\top x_3}^\infty
  \exp\{\sigma\left(\varrho_1 z_1+\varrho_3 z_3 \right)\}
  \phi\left(z_1,z_3;\rho_{13}\right)dz_1dz_3.
\end{array}
  \end{equation}

\item $E(y)$ for bivariate hurdle models \code{011} and \code{110} and
  univariate hurdle model \code{010} can be derived from that of
  trivariate model \code{111} by eliminating hurdles 1, 3, 1 and 3,
  respectively.  Likewise, the expectation of $y$ for bivariate hurdle
  models \code{100} and \code{001} can be derived from that of
  trivariate hurdle model \code{101} by eliminating hurdles 1 and 3,
  respectively.  Corresponding formulas of $E(y|y>0)= E(y)/P(y>0)$ for
  all this special cases implemented in \proglang{R} are presented in
  Table \ref{tab:typo}, using the following notations:
  $$
  \Psi_{2|1}=\rho_{12}\phi_1\Phi\left(\frac{\frac{\beta_2^\top x_2}{\sigma}
  -\rho_{12}\beta_1^\top x_1}{\sqrt{1-\rho_{12}^2}}\right)+\phi_2\Phi\left(\frac{\beta_1^\top x_1
  -\rho_{12}\frac{\beta_2^\top x_2}{\sigma}}{\sqrt{1-\rho_{12}^2}}\right),
  $$
  $$
  \Psi_{2|3}=\rho_{23}\phi_3\Phi\left(\frac{\frac{\beta_2^\top x_2}{\sigma}
  -\rho_{23}\beta_3^\top x_3}{\sqrt{1-\rho_{23}^2}}\right)+\phi_2\Phi\left(\frac{\beta_3^\top x_3
  -\rho_{23}\frac{\beta_2^\top x_2}{\sigma}}{\sqrt{1-\rho_{23}^2}}\right),
  $$
  where $\phi_1=\phi(\beta_1^\top x_1)$, $\phi_2\left(\frac{\beta_2^\top x_2}{\sigma}\right)$ and
  $\phi_3=\phi(\beta_3^\top x_3)$.

  Note that formulas of $E(y|y>0)$ for dependent trivariate hurdle models presented
  in Table \ref{tab:typo} are obtained by using closed forms of the
  following integrals :
  $$
  \int_{-\beta^\top x}^\infty\left[\rho z\Phi\left(\frac{\frac{\beta_2^\top x_2}{\sigma}+\rho z}
  {\sqrt{1-\rho^2}}\right)+\sqrt{1-\rho^2}\phi\left(\frac{\frac{\beta_2^\top x_2}{\sigma}+\rho z}
  {\sqrt{1-\rho^2}}\right)\right]\phi(z)dz
  $$
  $$
  =\rho\phi\left(\beta^\top x\right)\Phi\left(\frac{\frac{\beta_2^\top x_2}{\sigma}
  -\rho\beta^\top x}{\sqrt{1-\rho^2}}\right)+\phi\left(\frac{\beta_2^\top x_2}{\sigma}\right)
  \Phi\left(\frac{\beta^\top x-\rho\frac{\beta_2^\top x_2}{\sigma}}{\sqrt{1-\rho^2}}\right),
  $$
  $$
  \int_{-\beta^\top x}^\infty \exp\left\{\sigma\rho z\right\}\phi(z)dz=
  \exp\left\{\frac{\sigma^2\rho^2}{2}\right\}\Phi\left(\beta^\top x+\sigma\rho\right).
  $$

\end{itemize}

Denoting by $\hat{y}_i$ the fitted values of $y_i$ obtained by
estimating the best mean-square error predictor of $y_i$, namely
$E(y_i)$, with the maximum likelihood estimate of model parameters, we
define a pseudo coefficient of determination for a multiple hurdle
model using the following formula:

\begin{equation}\label{eq38}
R^2=1-\frac{RSS}{TSS},
\end{equation}

with $RSS=\sum (y_i - \hat{y}_i)^2$ the residual sum of squares and
$TSS=\sum (y_i - \hat{y}_0)^2$ the total sum of squares, where
$\hat{y}_0$ denotes the maximum likelihood estimate of $E(y_i)$ in the
multiple hurdle model without covariates (intercept-only model
\footnote{For multiple hurdle models involving many intercepts, the
estimation of a specification without covariates may face serious
numerical problems. If the \pkg{mhurdle} software fails to provide such an
estimate, the total sum of squares $TSS$ is computed by substituting
the sample average of $y$ for $\hat{y}_0$.}). Note that
this goodness of fit measure cannot exceed one but can be negative, as
a consequence of the non linearity of $E(y_i)$ with respect to the
parameters.

The extension of the McFadden likelihood ratio index for qualitative
response models to multiple hurdle models is straightforwardly
obtained by computing this index formula:

\begin{equation}\label{eq40}
\rho^2=1-\frac{\ln L(\hat{\theta})}{\ln L(\hat{\alpha})}=\frac{\ln
L(\hat{\alpha})-\ln L(\hat{\theta})}{\ln L(\hat{\alpha})},
\end{equation}

using the maximised log-likelihood functions of a multiple hurdle
model with covariates, $\ln L(\hat{\theta})$, and without covariates,
$\ln L(\hat{\alpha})$, respectively.  This goodness of fit measure
takes values within zero and one and, as it can be easily inferred
from the above second expression of $\rho^2$, it measures the relative
increase of the maximised log-likelihood function due to the use of
explanatory variables with respect to the maximised log-likelihood
function of a naive intercept-only model.

Model selection deals with the problem of discriminating between
alternative model specifications used to explain the same dependent
variable, with the purpose of finding the one best suited to explain
the sample of observations at hand.  This decision problem can be
tackled from the point of view of the model specification achieving
the best in-sample fit.

This selection criterion is easy to apply as it consists in comparing
one of the above defined measures of fit, computed for the competing
model specifications, after adjusting them for the loss of sample
degrees of freedom due to model parametrisation.  Indeed, the value of
these measures of fit can be improved by increasing model
parametrisation, in particular when the parameter estimates are
obtained by optimising a criteria functionally related to the selected
measure of fit, as is the case when using the $\rho^2$ fit measure
with a maximum likelihood estimate. Consequently, a penalty that
increases with the number of model parameters should be added to the
$R^2$ and $\rho^2$ fit measures to trade off goodness of fit
improvements with parameter parsimony losses.

To define an adjusted pseudo coefficient of determination, we rely on
\citet{THEIL/73}'s correction of $R^2$ in a linear regression model,
defined by

\begin{equation}\label{eq41}
\bar{R}^2=1-\frac{n-K_0}{n-K}\frac{RSS}{TSS},
\end{equation}

where $K$ and $K_0$ stand for the number of parameters of the multiple
hurdle model with covariates and without covariates, respectively
\footnote{When the mhurdle software fails to provide the parameter 
estimates of the intercept-only model and the total sum of squares 
$TSS$ is computed by substituting the sample average of $y$ for 
$\hat{y}_0$, $K_0$ is set equal to 1.}. Therefore, choosing the 
model specification with the largest $\bar{R}^2$ is equivalent to 
choosing the model specification with the smallest model residual 
variance estimate: $s^2=\frac{RSS}{n-K}$.

To define an adjusted likelihood ratio index, we replace in this
goodness of fit measure $\rho^2$ the log-likelihood criterion with the
Akaike information criterion $AIC=-2\ln
L(\hat{\theta})+2K$. Therefore, choosing the model specification with
the largest

\begin{equation}\label{eq42}
\bar{\rho}^2=1-\frac{\ln L(\hat{\theta})-K}{\ln L(\hat{\alpha})-K_0}
\end{equation}

is equivalent to choosing the model specification that minimises the
\citet{AKAIKE/73} predictor of the Kullback-Leibler Information
Criterion (KLIC). This criterion measures the distance between the
conditional density function $f(y|x;\theta)$ of a possibly
misspecified parametric model and that of the true unknown model,
denoted by $h(y|x)$. It is defined by the following formula:

\begin{equation}\label{eq43}
KLIC=E\left[\ln
  \left(\frac{h(y|x)}{f(y|x;\theta_\ast)}\right)\right]=\int \ln
\left(\frac{h(y|x)}{f(y|x;\theta_\ast)}\right)dH(y,x),
\end{equation}

where $H(y,x)$ denotes the distribution function of the true joint
distribution of $(y,x)$ and $\theta_\ast$ the probability limit, with
respect to $H(y,x)$, of $\hat{\theta}$ the so called quasi-maximum
likelihood estimator obtained by applying the maximum likelihood when
$f(y|x;\theta)$ is misspecified.


\subsection{Model selection using Vuong tests}

Model selection can also be tackled from the point of view of the model
specification that is favoured in a formal test comparing two model
alternatives.

This second model selection criterion relies on the use of a test
proposed by \citet{VUONG/89}. According to the rationale of this test,
the "best" parametric model specification among a collection of
competing specifications is the one that minimises the $KLIC$
criterion or, equivalently, the specification for which the quantity:

\begin{equation}\label{eq44}
E[\ln f(y|x;\theta_\ast)]=\int \ln f(y|x;\theta_\ast)dH(y,x)
\end{equation}

is the largest. Therefore, given two competing conditional models with
density functions $f(y|x;\theta)$ and $g(y|x;\pi)$ and parameter
vectors $\theta$ and $\pi$ of size $K$ and $L$, respectively, Vuong
suggests to discriminate between these models by testing the null
hypothesis:

$$
H_0 : E[\ln f(y|x;\theta_\ast)]=E[\ln g(y|x;\pi_\ast)]
\Longleftrightarrow
E\left[\ln\frac{f(y|x;\theta_\ast)}{g(y|x;\pi_\ast)}\right]=0,
$$

meaning that the two models are equivalent, against:

$$
H_f : E[\ln f(y|x;\theta_\ast)]>E[\ln g(y|x;\pi_\ast)]
\Longleftrightarrow
E\left[\ln\frac{f(y|x;\theta_\ast)}{g(y|x;\pi_\ast)}\right]>0,
$$

meaning that specification $f(y|x;\theta)$ is better than
$g(y|x;\pi)$, or against:

$$
H_g : E[\ln f(y|x;\theta_\ast)]<E[\ln g(y|x;\pi_\ast)]
\Longleftrightarrow
E\left[\ln\frac{f(y|x;\theta_\ast)}{g(y|x;\pi_\ast)}\right]<0,
$$

meaning that specification $g(y|x;\pi)$ is better than $f(y|x;\theta)$.

The quantity $E[\ln f(y|x;\theta_\ast)]$ is unknown but it can be
consistently estimated, under some regularity conditions, by $1/n$
times the log-likelihood evaluated at the quasi-maximum likelihood
estimator. Hence $1/n$ times the log-likelihood ratio (LR) statistic

\begin{equation}\label{eq45}
LR(\hat{\theta},\hat{\pi})=\sum_{i=1}^n
\ln\frac{f(y_i|x_i;\hat{\theta})}{g(y_i|x_i;\hat{\pi})}
\end{equation}

is a consistent estimator of
$E\left[\ln\frac{f(y|x;\theta_\ast)}{g(y|x;\pi_\ast)}\right]$. Therefore,
an obvious test of $H_0$ consists in verifying whether the LR
statistic differs from zero.  The distribution of this statistic can
be worked out even when the true model is unknown, as the
quasi-maximum likelihood estimators $\hat{\theta}$ and $\hat{\pi}$
converge in probability to the pseudo-true values $\theta_\ast$ and
$\pi_\ast$, respectively, and have asymptotic normal distributions
centred on these pseudo-true values.

The resulting distribution of $LR(\hat{\theta},\hat{\pi})$ depends on
the relation linking the two competing models. To this purpose, Vuong
differentiates among three types of competing models, namely: nested,
strictly non nested and overlapping.

A parametric model $G_\pi$ defined by the conditional density function
$g(y|x;\pi)$ is said to be nested in parametric model $F_\theta$ with
conditional density function $f(y|x;\theta)$, if and only if any
conditional density function of $G_\pi$ is equal to a conditional
density function of $F_\theta$ almost everywhere (disregarding any
zero probability sub-set of $(y,x)$ values, with respect to the true
distribution function $H(y,x)$). This means that we can write a
parametric constraint in the form $\theta=T(\pi)$, allowing to express
model $G_\pi$ as a particular case of model $F_\theta$. Within our
multiple hurdle special models this is the case when comparing two
specifications differing only with respect to the presence or the
absence of correlated disturbances.  For these models, it is
necessarily the case that $f(y|x;\theta_\ast)\equiv
g(y|x;\pi_\ast)$. Therefore $H_0$ is tested against $H_f$.

If model $F_\theta$ is misspecified, it has been shown by Vuong that:

\begin{itemize}
\item under $H_0$, the quantity $2LR(\hat{\theta},\hat{\pi})$ converges
  in distribution towards a weighted sum of $K+L$ iid $\chi^2(1)$
  random variables, where the weights are the $K+L$ almost surely real
  and non negative eigenvalues of the following $(K+L) \times (K+L)$
  matrix:

  $$
  W=\left[ \begin{array}{c}
  -B_f A_f^{-1}\quad \  \quad -B_{fg} A_g^{-1} \\
  B_{fg}^\top A_f^{-1}\quad \ \quad\quad B_g A_g^{-1}
  \end{array} \right],
  $$

  where

  $$
  A_f=E\left(\frac{\partial^2\ln
      f(y|x;\theta_\ast)}{\partial\theta\partial\theta^\top}\right),\quad
  A_g=E\left(\frac{\partial^2\ln
      g(y|x;\pi_\ast)}{\partial\pi\partial\pi^\top}\right),
  $$
  $$
  B_f=E\left(\frac{\partial\ln
      f(y|x;\theta_\ast)}{\partial\theta}\frac{\partial\ln
      f(y|x;\theta_\ast)}{\partial\theta^\top}\right),\quad
  B_g=E\left(\frac{\partial\ln
      g(y|x;\pi_\ast)}{\partial\pi}\frac{\partial\ln
      g(y|x;\pi_\ast)}{\partial\pi^\top}\right),
  $$
  $$
  B_{fg}=E\left(\frac{\partial\ln
      f(y|x;\theta_\ast)}{\partial\theta}\frac{\partial\ln
      g(y|x;\pi_\ast)}{\partial\pi^\top}\right).
  $$

  To simplify the computation of this limiting distribution,
  one can alternatively use the weighted sum of $K$ iid
  $\chi^2(1)$ random variables, where the weights are the $K$
  almost surely real and non negative eigenvalues of the following
  smaller $K \times K$ matrix:

  $$
  \underline{W}=B_f\left[DA_g^{-1}D^\top -A_f^{-1}\right],
  $$

  where $D=\frac{\partial T(\pi_\ast)}{\partial\pi^\top}$.

\item under $H_f$, the same statistic converge almost surely towards
  $+\infty$.

\end{itemize}

Performing this standard LR test for nested models, requires to
replace the theoretical matrices $W$ and $\underline{W}$ by a
consistent estimator. Such an estimator is obtained by substituting
matrices $A_f$, $A_g$, $B_f$, $B_g$ and $B_{fg}$ for their sample
analogue:

$$
\hat{A}_f=\frac{1}{n}\sum_{i=1}^n\frac{\partial^2\ln
  f(y_i|x_i;\hat{\theta})}{\partial\theta\partial\theta^\top},\quad
\hat{A}_g=\frac{1}{n}\sum_{i=1}^n\frac{\partial^2\ln
  g(y_i|x_i;\hat{\pi})}{\partial\pi\partial\pi^\top},
$$
$$
\hat{B}_f=\frac{1}{n}\sum_{i=1}^n\frac{\partial\ln
  f(y_i|x_i;\hat{\theta})}{\partial\theta}\frac{\partial\ln f(y_i|x_i;
  \hat{\theta})}{\partial\theta^\top},\quad
\hat{B}_g=\frac{1}{n}\sum_{i=1}^n\frac{\partial\ln
  g(y_i|x_i;\hat{\theta})}{\partial\theta}\frac{\partial\ln g(y_i|x_i;
  \hat{\theta})}{\partial\theta^\top},
$$
$$
\hat{B}_{fg}=\frac{1}{n}\sum_{i=1}^n\frac{\partial\ln
  f(y_i|x_i;\hat{\theta})}{\partial\theta}\frac{\partial\ln g(y_i|x_i;
  \hat{\theta})}{\partial\theta^\top}
$$

and $D$ for $\hat{D}=\partial T(\hat{\pi})/\partial\pi^\top$.

The density function of this asymptotic test statistic has not
been worked out analytically. Therefore, we compute it by simulation.

Hence, for a test with critical value $c$, $H_0$ is rejected in favour
of $H_f$ if $2LR(\hat{\theta},\hat{\pi})>c$ or if the p-value
associated to the observed value of $2LR(\hat{\theta},\hat{\pi})$ is
less than the significance level of the test.

Note that, if model $F_\theta$ is correctly specified, the
asymptotic distribution of the LR statistic is, as expected, a
$\chi^2$ random variable with $K-L$ degrees of freedom.

Two parametric models $F_\theta$ and $G_\pi$ defined by conditional
distribution functions $f(y|x;\theta)$ and $g(y|x;\pi)$ are said to be
strictly non-nested, if and only if no conditional distribution
function of model $F_\theta$ is equal to a conditional distribution
function of $G_\pi$ almost everywhere, and conversely. Within multiple
hurdle special models this is the case when comparing two
specifications differing with respect either to the censoring
mechanisms in effect or to the functional form of the desired
consumption equation. For these models, it is necessarily the case
that $f(y|x;\theta_\ast)\neq g(y|x;\pi_\ast)$ implying that both
models are misspecified under $H_0$.

For such strictly non-nested models, Vuong has shown that:

\begin{itemize}
\item under $H_0$, the quantity $n^{-1/2}LR(\hat{\theta},\hat{\pi})$
  converges in distribution towards a normal random variable with zero
  expectation and variance:
 $$
 \omega^2=V\left(\ln\frac{f(y|x;\theta_\ast)}{g(y|x;\pi_\ast)}\right)
 $$
 computed with respect to the distribution function of the true joint
 distribution of $(y,x)$.
\item under $H_f$, the same statistic converge almost surely towards
  $+\infty$.
\item under $H_g$, the same statistic converge almost surely towards
  $-\infty$.
\end{itemize}

Hence, $H_0$ is tested against $H_f$ or $H_g$ using the standardised
LR statistic:

\begin{equation}\label{eq46}
T_{LR}=\frac{LR(\hat{\theta},\hat{\pi})}{\sqrt{n}\hat{\omega}},
\end{equation}

where $\hat{\omega}^2$ denotes the following strongly consistent
estimator for $\omega^2$:

$$
\hat{\omega}^2=\frac{1}{n}\sum_{i=1}^n
\left(\ln\frac{f(y_i|x_i;\hat{\theta})}{g(y_i|x_i;\hat{\pi})}\right)^2
-\left(\frac{1}{n}\sum_{i=1}^n
  \ln\frac{f(y_i|x_i;\hat{\theta})}{g(y_i|x_i;\hat{\pi})}\right)^2 .
$$

As a consequence, for a test with critical value $c$, $H_0$ is
rejected in favour of $H_f$ if $T_{LR}>c$ or if the p-value associated
to the observed value of $T_{LR}$ in less than the significance level
of the test. Conversely, $H_0$ is rejected in favour of $H_g$ if
$T_{LR}<-c$ or if the p-value associated to the observed value of
$|T_{LR}|$ in less than the significance level of the test.

Note that, if one of models $F_\theta$ or $G_\pi$ is assumed to be
correctly specified, the \citet{COX/61,COX/62} LR test of non nested
models needs to be used. Because this test is computationally awkward
to implement and not really one of model selection, as it can lead to
reject both competing models, it has not been programmed in
\pkg{mhurdle}.

Two parametric models $F_\theta$ and $G_\pi$ defined by conditional
distribution functions $f(y|x;\theta)$ and $g(y|x;\pi)$ are said to be
overlapping, if and only if part of the conditional distribution
function of model $F_\theta$ is equal to the conditional distribution
function of $G_\pi$ but none of these models is nested in the
other. Within multiple hurdle special models this is the case when
comparing two specifications differing only with respect to the
covariates taken into consideration, some of them being common to
both models and others specific. For these models it is not clear
\emph{a priori} as to whether or not $f(y|x;\theta_\ast)=
g(y|x;\pi_\ast)$ almost everywhere, except if we know \emph{a priori}
that at least one of the two competing models is correctly
specified. As a consequence, the form of the asymptotic distribution
of $LR(\hat{\theta},\hat{\pi})$ under $H_0$ is unknown, which prevents
from performing a model selection test based on this statistic.

In the general case where both competing models are wrongly specified,
Vuong suggests a sequential procedure which consists in testing first
whether or not the variance $\omega^{2}$ equals zero (since
$f(y|x;\theta_\ast)= g(y|x;\pi_\ast)$ almost everywhere if and only if
$\omega^{2}=0$) and then, according to the outcome of this test, in
using the appropriate asymptotic $LR(\hat{\theta},\hat{\pi})$
distribution to perform the model selection test.

To test $H_0^{\omega}: \omega^2=0$ against $H_A^{\omega}: \omega^2\neq
0$, Vuong suggests to use, as a test statistic, the above defined
strongly consistent estimator for $\omega^2$, $\hat{\omega}^2$, and
proves that:

\begin{itemize}
\item under $H_0^{\omega}$, the quantity $n\hat{\omega}^2$ converges
  in distribution towards the same limiting distribution like that
  of statistic $2LR(\hat{\theta},\hat{\pi})$ when used for discriminating
  two misspecified nested models.

\item under $H_A^{\omega}$, the same statistic converge almost surely towards
  $+\infty$.
\end{itemize}

Therefore, performing this variance test requires to compute the
eigenvalues of a consistent estimate of matrix $W$ or $\underline{W}$,
and derive by simulation the density function of the corresponding
weighted sum of iid $\chi^2(1)$ random variables.

Hence, for a test with critical value $c$, $H_0^\omega$ is rejected in
favour of $H_A^\omega$ if $n\hat{\omega}^2>c$ or if the p-value
associated to the observed value of $n\hat{\omega}^2$ is less than the
significance level of the test.

Note, that an asymptotically equivalent test is obtained by replacing
in statistics $n\hat{\omega}^2$, $\hat{\omega}^2$ by:
$$
\tilde{\omega}^2=\frac{1}{n}\sum_{i=1}^n
\left(\ln\frac{f(y_i|x_i;\hat{\theta})}{g(y_i|x_i;\hat{\pi})}\right)^2.
$$

The second step in discriminating two overlapping models depends on
the outcome of the variance test.

\begin{itemize}
\item If $H_0^\omega$ is not rejected, one should conclude that the
  two models cannot be discriminated given the data, since assuming
  $\omega^2=0$ implies that $H_0$ means that the two models are
  equivalent.
\item If $H_0^\omega$ is rejected, the test of $H_0$ against $H_f$ or
  $H_g$ must be carried out using the standardised LR statistic
  $T_{LR}$, as for discriminating between two strictly non-nested
  models. Indeed, $H_0$ is still possible when $\omega^2\neq0$. Note,
  that this sequential procedure of testing $H_0$ against $H_f$ or
  $H_g$ has a significance level bounded above by the maximum of the
  significance levels used for performing the variance and the
  standardised LR tests.
\end{itemize}

Finally, if one of the two competing models is supposed to be
correctly specified, then the two models are equivalent if and only
if the other model is correctly specified and if and only if the
conditional density functions of the two models are identical almost
everywhere. In this case we can bypass the variance test and directly
construct a model selection test based on the
$2LR(\hat{\theta},\hat{\pi})$ test statistic used for discriminating
between two nested models.


\section{Software rationale}
\label{sec:software}

There are three important issues to be addressed to correctly
implement in \proglang{R} the modelling strategy described in the
previous sections. The first one is to provide a good interface to
describe the model to be estimated. The second one is to find
good starting values for computing model estimates. The third
one is to have flexible optimisation tools for likelihood
maximisation.

\subsection{Model syntax}
\label{sec:modeldesc}


In \proglang{R}, the model to be estimated is usually described using
formula objects, the left-hand side denoting the censored dependent
variable \texttt{y} and the right-hand side the functional relation
explaining \texttt{y} as a function of covariates. For example,
\texttt{y \~{} x1 + x2 * x3} indicates that \texttt{y} linearly depends
on variables \texttt{x1}, \texttt{x2}, \texttt{x3} and on the
interaction term \texttt{x2} times \texttt{x3}.

For the models implemented in \pkg{mhurdle}, four kinds of covariates
should be specified: those of 
\begin{itemize}
\item the good selection equation (hurdle 1) denoted $x_1$,
\item the desired consumption equation (hurdle 2), denoted $x_2$,
\item the purchasing equation (hurdle 3), denoted $x_3$,
\item the variance equation, denoted $x_4$.
\end{itemize}

To define a model with several kinds of covariates, a general solution
is given by the \pkg{Formula} package developed by
\citet{ZEIL/CROI/10}, which provides extended formula objects. To
define a model where \texttt{y} is the censored dependent variable,
\texttt{x11} and \texttt{x12} two covariates for the good selection
equation, \texttt{x21} and \texttt{x22} two covariates for the desired
consumption equation, and \texttt{x31} and \texttt{x32} two covariates
for the purchasing equation, we use the following commands :

<<label = Formula>>=
library("Formula")
f <- Formula(y ~ x11 + x12  | x21 + x22 | x31 + x32)
@

\subsection{Starting values and optimisation}
\label{sec:startvalues}

For the models we consider, the log-likelihood function will be, in
general, not concave. Moreover, this kind of models are highly non
linear with respect to parameters, and therefore difficult to
estimate. For these reasons, the question of finding good starting
values for the iterative computation of parameter estimates is
crucial.

As a less computer intensive alternative to maximum likelihood
estimation, \citet{HECKM/76} has suggested a two step estimation
procedure based on a respecification of the censored variable linear
regression model, sometimes called ``Heckit'' model, avoiding
inconsistency of the ordinary least-squares estimator.  This two step
estimator is consistent but inefficient. It is implemented in package
\pkg{sampleSelection} \citep{TOOME/HENNI/08}.

According to \citet{CARL/CROI/HOAR/08} experience in applying this
estimation procedure to double hurdle models, this approach doesn't
seem to work well with correlated hurdle models. Indeed, except for
the very special case of models \code{100}, \code{010} and \code{001},
the probability of observing a censored purchase is not that of a
simple probit model (see Table \ref{tab:typo}).

As noted previously, for uncorrelated single hurdle models, the
estimation may be performed in a sequence of two simple estimations,
namely the maximum likelihood estimation of a standard dichotomous
probit model, followed by the ordinary least-squares estimation of a
linear, log-linear or linear-truncated regression model. In the last
case, package \pkg{truncreg} \citep{TRUNCR/09} is used.

%% For correlated single hurdle 1 model \code{100}, the maximum
%% likelihood estimate of the parameters of the corresponding
%% uncorrelated model ($\rho_{12}=0$) is used as starting values.

%% For P-Tobit models (\code{001} and \code{011}), the starting values
%% are computed using an Heckman-like two step procedure. In the first
%% step, parameters $\beta_3$ are estimated using a simple probit. In the
%% second step, a linear regression model is estimated by ordinary least
%% squares using the sub-sample of uncensored observations and, as
%% dependent variable: $y_i\Phi(\hat{\beta}_3^{\top}x_{3i})$ in the case
%% of a normal specification and $\ln y_i +\ln
%% \Phi(\hat{\beta}_3^{\top}x_{3i})$ in the case of a log-normal
%% specification.

%% For Tobit model \code{010}, the least squares estimate of the linear regression
%% model is used as starting values.

%% For double hurdle model \code{110}, the starting values for $\beta_1$
%% are obtained by estimating a probit model and those for $\beta_2$
%% using a least squares estimate with the truncated sample of a linear
%% regression model assuming $\rho_{12}=0$.

%% Finally, for models involving hurdles 1 and 3 (models \code{101} and
%% \code{111}), we use two probit models to get starting values for
%% $\beta_1$ and $\beta_2$. Then, we estimate a linear regression model
%% by ordinary least squares with the sub-sample of uncensored
%% observations, as dependent variable, and assuming no correlation
%% between the desired consumption equation and these two hurdles:
%% $y_i\Phi(\hat{\beta}_3^{\top}x_{3i})$, in the case of a normal
%% specification and $\ln y_i +\ln \Phi(\hat{\beta}_3^{\top}x_{3i})$ in
%% the case of a log-normal specification.

%% \subsection{Optimisation}
%% \label{sec:optimisation}

Two kinds of algorithms are currently used for maximum likelihood
estimation. The first kind of algorithms can be called ``Newton-like''
methods.  With these algorithms, at each iteration, the Hessian matrix
of the log-likelihood is computed, using either the second derivatives
of the log-likelihood (Newton-Raphson method) or the outer product of
the gradient (Berndt, Hall, Hall, Hausman or \textsc{bhhh}
method). This approach is very powerful if the log-likelihood is
well-behaved, but it may perform poorly otherwise and fail after a few
iterations.

The second algorithm, called Broyden, Fletcher, Goldfarb, Shanno or
\textsc{bfgs} method, updates at each iteration an estimate of the
Hessian matrix of the log-likelihood. It is often more robust and may
perform better in cases where the formers don't work.

Two optimisation functions are included in core \proglang{R}:
\code{nlm}, which uses the Newton-Raphson method, and \code{optim} ,
which uses the \textsc{bfgs} method (among others). The recently
developed \pkg{maxLik} package by \citet{MAXLIK/08} provides a unified
framework. With a unified interface, all the previously described
methods are available.

The behaviour of \code{maxLik} can be controlled by the user using
\code{mhurdle} arguments like \texttt{print.level} (from 0-silent to
2-verbal), \texttt{iterlim} (the maximum number of iterations),
\texttt{methods} (the method used, one of \texttt{"nr"},
\texttt{"bhhh"} or \texttt{"bfgs"}) that are passed to \code{maxLik}.

Some models require the computation of the bivariate normal cumulative
density function. We use the \pkg{pbivnorm} package \citep{KENKE:11}
which provides a vectorised (and therefore fast and convenient)
function to compute the bivariate normal cdf.


\section{Examples}
\label{sec:examples}

The package is loaded using:

<<label = Library>>=
library("mhurdle")
@

To illustrate the use of \pkg{mhurdle}, we use one surveys conducted
by the Bureau of Labour Statistics of the U.S. Department of Labour,
called the "Interview Survey". Data from 25813 households on all
expenditures are collected, on a quarterly basis. The micro-data files
are publicly available on the website of the Bureau of Labour
Statistics, and may be downloaded and used without permission.  We use
a small subset of 1000 randomly selected households. The total
expenditure is divided in 14 main chapters (food, alcool, housing,
apparel, transport, health, entertainment, personal care, reading,
eduction, tobacco, miscellanous expenditures, cash contributions and
insurance) and we also report expenditures on fees and admissions
(called \code{shows}) which is part of the entertainment chapter and
trip expenditures (called \code{vacations} on the data set).

All the expenditures are in tousands of USD, measured on an annual
basis and are divided by the number of consumption units (obtained by
counting for one the first adult of the household, 0.7 the subsequent
adults and 0.5 every other person aged under 18).


<<eval = TRUE, label = Data>>=
data("Interview", package = "mhurdle")
head(Interview, 3)
mean(Interview$shows == 0)
max(Interview$shows)
@

The covariates are :

\begin{description}
\item[income] the anual net income by consumption unit,
\item[smsa]{does the household live in a SMSA (\code{yes} or
    \code{no}),}
\item[age]{the age of the reference person of the household,}
\item[educ]{the number of year of education of the reference person of
    the household,}
\item[sex]{the sex of the reference person of the household
    (\code{male} and \code{female}),}
\item[size]{the number of persons in the household,}
\item[month]{the month of the interview (between 1 and 12),}
\end{description}


\subsection{Estimation}

The estimation is performed using the \code{mhurdle} function, which
has the following arguments:

\begin{description}
\item[formula:] a formula describing the model to estimate. It should
  have between two and four parts on the right-hand side specifying,
  in the first part, the good selection equation covariates, in the
  second part, the desired consumption equation covariates, in the
  third part, the purchasing equation covariates and in the fourth
  part, the covariates of the variance equation.
\item[data:] a data frame containing the observations of the variables
  present in the formula.
\item[subset, weights, na.action:] these are arguments passed on to
  the \code{model.frame} function in order to extract the data
  suitable for the model. These arguments are present in the \code{lm}
  function and in most of the estimation functions.
\item[start:] the starting values. If \code{NULL}, the starting values
  are computed as described in section 4.2.
\item[dist:] this argument indicates the functional form of the
  desired consumption equation, which may be either log-normal
  \code{"ln"} (the default), two-parameters log-normal \code{"ln2"},
  normal \code{"n"}, truncated normal \code{"tn"}, Box-Cox
  \code{"bc"}, two-parameters Box-Cox \code{"bc2"} or inverse
  Hyperbolic Sine \code{"ihs"},
\item[scaled:] if \code{TRUE}, the dependent variable is divided by
  the geometric mean of the positive values,
\item[corr:] this boolean argument indicates whether the disturbance
  of the different equations are correlated, the default value is
  \code{FALSE},
\item[robust:] if \code{TRUE}, transformations of some parameters are
  used, so that they lie in the required range (positive values for
  the standard deviation and for the position parameter, between -1
  and +1 for the coefficients of correlation),
\item[...] further arguments that are passed to the optimisation
  function \code{maxLik}.
\end{description}

\subsubsection{One equation models}

We start with models that only contain the consumption equation. In
this case, the only source of null consumption is the lack of
ressources. In this case, the distribution of $y^*$ must admit
negative values, this is the case for the two-parameters box-cox
normal distribution, the log-normal distribution with a position
parameter and the normal distribution (this later case corresponds to
the standard tobit model). We use the expenditures on vacation as an
exemple and the results are presented in table~\ref{tab:oneq}, using
the \pkg{texreg} library.

<<label = onequation>>=
Sn <- mhurdle(vacations ~ 0 | I(car > 0) + size + linc + linc2 + age + age2,
              Interview, dist = "n", h2 = TRUE, corr = FALSE,
              method = "bhhh", print.level = 0)
Sbc <- update(Sn, dist = "bc")
Sln <- update(Sn, dist = "ln")
library(texreg)
@ 

<<eval = FALSE, echo = TRUE, results = hide>>=
texreg(list(Sn, Sbc, Sln), 
       custom.model.names = c("normal tobit", "log-normal tobit", "box-cox tobit"), 
       caption = "Estimation of one-equation models", 
       label = "tab:oneq", pos = "ht")
@ 

<<echo = FALSE, results = tex>>=
texreg(list(Sn, Sbc, Sln), 
       custom.model.names = c("normal tobit", "log-normal tobit", "box-cox tobit"), 
       caption = "Estimation of one-equation models", 
       label = "tab:oneq", pos = "ht")
@ 

\subsubsection{Simple selection models}

The original independent single hurdle models proposed by
\citet{CRAGG/71} required the distinct estimation of two models : a
probit to explain the null \emph{vs} the positive observations and an
estimation on the truncated sample in order to explain the level of
consumption of those who consume the good. This latter estimation can
be realised with the hypothesis that the distribution of the errors
are truncated normal, log-normal or box-cox normal. We use the
expenditure on food-away as an exemple and the results are presented
in table~\ref{tab:sep}.

<<label = selection>>=
Stn <- mhurdle(foodaway ~ size + smsa + age + age2 | linc + linc2, Interview, 
               dist = "n", h2 = FALSE, corr = FALSE, method = "bhhh", print.level = 0)
Sbc <- update(Stn, dist = "bc")
Sln <- update(Stn, dist = "ln")
@ 

%% <<label = selection>>=
%% Sbc <- mhurdle(foodaway ~ size + smsa + age + age2 | linc + linc2, Interview, 
%%                dist = "bc", h2 = FALSE, corr = FALSE, method = "bhhh", print.level = 0)
%% Sln <- update(Sbc, dist = "ln")
%% @ 


The dependent model is easily obtained by setting the \code{corr}
argument to \code{TRUE}. For the log-normal model, we obtain :


<<>>=
Slnd <- update(Sln, corr = TRUE)
coef(summary(Slnd), "corr")
@ 


The coefficient of correlation between the selection and the
consumption is in this case highly significant.


<<eval = FALSE, echo = TRUE, results = hide>>=
texreg(list(Stn, Sln, Sbc), 
       custom.model.names = c("truncated-normal", "log-normal", "box-cox"), 
       caption = "Estimation of single hurdle selection models", 
       label = "tab:sep", pos = "ht")
@ 


<<echo = FALSE, results = tex>>=
texreg(list(Stn, Sln, Sbc), 
       custom.model.names = c("truncated-normal", "log-normal", "box-cox"), 
       caption = "Estimation of single hurdle selection models", 
       label = "tab:sep", pos = "ht")
@ 



\subsubsection{P-tobit model}

Apparel is a good candidate to illustrate the estimation of a P-tobit
model as the infrequency of purchase is clearly the only source of
null expenditure in this case. To explain the probability of
purchasing during the quarter of the survey, we use the month (as
expenditures in apparels can be concentrated on certain periods during
the year because of the sales) and the fact that the household lives
in a \textsc{smsa}. No corner solution is allowed as the good is
consumed even for households with low ressources, so we choose a
log-normal distribution.

<<ptobit>>=
ptobit <- mhurdle(apparel ~ 0 | linc + linc2 | factor(month) + smsa,
                  Interview, corr = TRUE, dist = "ln", h2 = FALSE,
                  method = "bhhh")
@ 



To illustrate the use of \code{mhurdle}, we estimate several models
explaining the expenditure on fees and admission. The expenditure is
supposed to depends on income and its square, age and its square,
education, the size of the household and the fact that the household
lives in a smsa. We first estimate a triple hurdle model, using
\code{educ} and \code{size} as covariates for the selection equation
and \code{age} and \code{smsa} as covariates for the purchasing
equation. We use a two-parameter log-normal distribution and a general
structure of correlation between the errors of the three equations is
estimated.

<<label = ThreeHurdlesDep>>=
H3D <- mhurdle(shows ~ educ + size | linc + linc2 | age + age2 + smsa, 
               Interview, dist = "ln", h2 = TRUE, corr = TRUE, method = "bhhh")
@ 

The independent version of the same model is obtained by setting the
\code{corr} argument to \code{FALSE}.

<<label = ThreeHurdlesInd>>=
H3I <- update(H3D, corr = FALSE)
@ 

A three equation - two hurdles model is obtained by making the second
hurdle ineffective,

<<label = TwoHurdlesThreeEqs>>=
H2D <- update(H3D, h2 = FALSE)
@ 

We then estimate a double-hurdle selection model by using a two-part
formula with \code{educ}, \code{size}, \code{age} and \code{smsa} in
the first part which describe the selection process:

<<label = Selection>>=
S2D <- mhurdle(shows ~ educ + size + age + age2 + smsa | linc + linc2, 
               Interview, dist = "ln", h2 = TRUE, corr = TRUE, method = "bhhh")
@ 

Finaly, we estimate a double-hurdle p-tobit model, all the covariates
which were in the first part of the formula being moved to the third
part:

<<label = Ptobit>>=
P2D <- mhurdle(shows ~ 0 | linc + linc2 | educ + size + age + age2 + smsa, 
               Interview, dist = "ln", h2 = TRUE, corr = TRUE, method = "bhhh")
@ 

\subsection{Methods}

A \code{summary} method is provided for \code{mhurdle} objects :


<<eval = TRUE, label = Summary>>=
summary(H3D)
@

This method displays the percentage of 0 in the sample, the
table of parameter estimates, and two measures of goodness of fit.

\code{coef}, \code{vcov}, \code{logLik}, \code{fitted} and
\code{predict} methods are provided in order to extract part of the
results.

Parameter estimates and the estimated asymptotic variance matrix of maximum
likelihood estimators are extracted using the usual \code{coef} and
\code{vcov} functions. \code{mhurdle} object methods have a second
argument indicating which subset has to be returned (the default is to
return all).

<<eval = TRUE, label = Methods>>=
coef(H3D, "h2")
coef(H3D, "h1")
coef(H3D, "sd")
coef(summary(H3D), "h3")
vcov(H3D, "h3")
@

Log-likelihood may be obtained for the estimated model or for a
``naive'' model, defined as a model without covariates :


<<eval = TRUE, label = LogLik>>=
logLik(H3D)
logLik(H3D, naive = TRUE)
@

Fitted values are obtained using the \code{fitted} method.  The output
is a matrix whose two columns are the estimated probability of
censoring $\mbox{P}(y=0)$ and the estimated expected value of an
uncensored dependent variable observation $\mbox{E}(y|y>0)$.


<<eval = FALSE, label = Fitted>>=
head(fitted(H3D))
@

A \code{predict} function is also provided, which returns the same
two columns for given values of the covariates.


<<eval = FALSE, label = Predict>>=
predict(H3D,
        newdata = data.frame(
            comics = c(0, 1, 2),
            gender = c("female", "female", "male"),
            age = c(20, 18, 32),
            educ = c(10, 20, 5),
            incum = c(4, 8, 2),
            size = c(2, 1, 3)))
@

For model evaluation and selection purposes, goodness of fit measures
and Vuong tests described in section 3 are provided. These criteria
allow to select the most empirically relevant model specification.

Two goodness of fit measures are provided. The first measure is an
extension to limited dependent variable models of the classical
coefficient of determination for linear regression models. This pseudo
coefficient of determination is computed both without (see formula
\ref{eq38}) and with (see formula \ref{eq41}) adjustment for the loss
of sample degrees of freedom due to model parametrisation.  The
unadjusted coefficient of determination allows to compare the goodness
of fit of model specifications having the same number of parameters,
whereas the adjusted version of this coefficient is suited for
comparing model specifications with a different number of parameters.


<<eval = FALSE, label = Rsq1>>=
rsq(H3D, type = "coefdet")
@

The second measure is an extension to limited dependent variable
models of the likelihood ratio index for qualitative response
models. This pseudo coefficient of determination is also computed both
without (see formula \ref{eq40}) and with (see formula \ref{eq42})
adjustment for the loss of sample degrees of freedom due to model
parametrisation, in order to allow model comparisons with the same
or with a different number of parameters.


<<eval = TRUE, label = Rsq2>>=
rsq(H3D, type = "lratio", adj = TRUE)
@

The Vuong test based on the $T_{LR}$ statistic, as presented in
section 3.3 (see formula \ref{eq46}), is also provided as a criteria
for model selection within the family of 8 strictly non-nested models
of Figure 1\footnote{Note that Vuong tests for strictly non-nested
  models can be performed using the \code{vuong} function of the
  \pkg{pscl} package of \citet{JACKM:11} for \code{glm} models and
  some specific count data models.} :


<<eval = TRUE, label = VuongUnnested>>=
vuongtest(S2D, P2D)
@

According to this outcome, the null hypothesis stating the equivalence
between the two models is rejected in favour of the alternative
hypothesis stating that \code{P2D} is better than \code{S2D}.

Testing the hypothesis of no correlation between the good selection
mechanism, the purchasing mechanism and the desired consumption
equation can be performed as a Vuong test of selection between two
nested models, differing only with respect to the value of the
correlation coefficients, namely the test of the hypothesis
$H_0: \rho_{12}=\rho_{13}=\rho_{23}=0$, specifying an independent
mhurdle model, against the alternative hypothesis specifying a
corresponding dependent mhurdle model. This test is performed using
the log-likelihood ratio (LR) statistic (\ref{eq45}). As explained in
section 3.3, the critical value or the p-value to be used to perform
this test is not the same depending on the model builder believes or
not that his unrestricted model is correctly specified.  In the first
case, the p-value is computed using the standard chi square
distribution, whereas in the second case a weighted chi square
distribution is used.

<<eval = TRUE, label = VuongNested>>=
vuongtest(H3D, H3I, type = 'nested', hyp = TRUE)
vuongtest(H3D, H3I, type = 'nested', hyp = FALSE)
@

According to these outcomes, the null hypothesis of zero correlation
is rejected if the unrestricted model is assumed to be correctly
specified and accepted otherwise.

%% ## Testing this hypothesis of no correlation by assuming the unrestricted
%% ## model correctly specified, can be also performed by means of the
%% ## classical Wald test, using the t-statistic or the p-value of the
%% ## correlation coefficient estimate presented in the table of parameter
%% ## estimates of the dependent model (\code{m100d}).

%% <<eval = TRUE>>=
%% coef(summary(m100d), "corr")
%% @

%% According to this test outcome, the hypothesis of zero correlation
%% is accepted at a little less stringent significance level than with
%% a Vuong test.

Finally, to illustrate the use of the Vuong test for discriminating
between two overlapping models, we consider a slightly different
selection model obtained by removing the \code{age} covariate and
adding the \code{sex} covariate :

<<eval = TRUE, label = SelectionBis>>=
S2Db <- mhurdle(shows ~ educ + size + sex + smsa | linc + linc2, 
                Interview, dist = "ln", h2 = TRUE, corr = TRUE, method = "bhhh")
@ 

In this case, the Vuong test is performed in two steps. Firstly a test
of the null hypothesis $\omega^2=0$, meaning that the two models are
equivalent, is undertaken.

<<eval = TRUE, label = VuongOverlap>>=
vuongtest(S2D, S2Db, type="overlapping")
@

The null hypothesis is not rejected. 

%% This null hypothesis is here strongly rejected. Therefore, we can test
%% the equivalence of these two models as if they were strictly
%% non-nested.

%% <<eval = TRUE>>=
%% vuongtest(m010, m010bis, type="non-nested")
%% @

%% According to the outcome of this second test, we conclude that these two
%% model specifications cannot be empirically discriminated.

If one of two overlapping models is assumed to be correctly specified,
we can bypass the first step of this Voung test (the variance test) and
proceed as if we had to discriminate between two nested models.

<<eval = TRUE, label = VuongOverlap2>>=
vuongtest(S2D, S2Db, type="overlapping", hyp=TRUE)
@

%Once again, the equivalence of the two models is not rejected.

\section{Conclusion}
\label{sec:conclusion}

\pkg{mhurdle} aims at providing a unified framework allowing to
estimate and assess a variety of extensions of the standard Tobit
model particularly suitable for single-equation demand analysis not
currently implemented in \proglang{R}. It explains the presence of a
large proportion of zero observations for a dependent variable by
means of up to three censoring mechanisms, called hurdles.  Inspired
by the paradigms used for analysing censored household expenditure
data, these hurdles express: (i) a non economic decision mechanism for
a good rejection or selection motivated by ethical, psychological or
social considerations; (ii) an economic decision mechanism for the
desired level of consumption of a previously selected good, which can
turn out to be negative leading to a nil consumption; (iii) an
economic or non economic decision mechanism for the time frequency at
which the desired quantity of a selected good is bought or
consumed. Unexplained interdependence between latent variables is
modelled by assuming a possible correlation between the random
disturbances in the model relations. Despite the particular area of
application from which the above mentioned censoring mechanisms stem,
the practical scope of \pkg{mhurdle} models doesn't seem to be
restricted to empirical demand analysis.

To provide an operational and efficient statistical framework,
\pkg{mhurdle} models are specified in a fully parametric form allowing
statistical estimation and testing within the maximum likelihood
inferential framework. Tools for model evaluation and selection are
provided, based on the use of goodness of fit measure extensions of
the classical coefficient of determination and of the likelihood ratio
index of McFadden, as well as on the use of Vuong tests for nested,
strictly non-nested and overlapping model comparison when none, one or
both of two competing models are misspecified.

Tests of \pkg{mhurdle} computing procedures with a wide variety of
simulated and observational data have proved the performance and
robustness of \pkg{mhurdle} package. Still, extensions and
improvements of the software are \change{continuing}.

%% Other desirable extensions, like the use of more general functional
%% forms of the desired consumption relation\footnote{See
%%   \citet{POIR/78}, \citet{LANK/WYCK/91} and \citet{JONE/YEN/00}.} or
%% of less stringent distributional assumptions on which semi-parametric
%% or nonparametric estimation methods are based, will be tackled once
%% the actual scope of our models is established through diversified
%% empirical applications. Research is continuing in this direction.



\bibliography{bibliomhurdle}

\end{document}



  \begin{tabular}{|l|lll|lll|l|l|l|}\hline\hline
id & $h_1$ & $h_2$ & $h_3$ & $\rho_{12}$ & $\rho_{13}$ & $\rho_{23}$  & $P(y=0)$ & $f_+(y)$ & $\mbox{E}(y \mid y > 0)$\\\hline\hline
1 &$\square$&$\square$&$\square$&$\square$&$\square$&$\square$& $0$ & $\frac{1}{\sigma y}\phi\left(\frac{\ln y - \beta_2^{\top}x_2}{\sigma}\right)$ & $\exp\left\{\beta_2^{\top}x_2 + \frac{\sigma^2}{2}\right\}$\\
2i &$\blacksquare$&$\square$&$\square$&$\square$&$\square$&$\square$& $1 - \Phi_1$& $\frac{1}{\sigma y}\phi\left(\frac{\ln y -\beta_2^{\top}x_2}{\sigma}\right)\Phi_1$& $\exp\left\{\beta_2^{\top} x_2 + \frac{\sigma^2}{2}\right\}$\\
2d &$\blacksquare$&$\square$&$\square$&$\blacksquare$&$\square$&$\square$ & $1 - \Phi_1$& $\frac{1}{\sigma y}\phi\left(\frac{\ln y -\beta_2^{\top}x_2}{\sigma}\right)\Phi\left(\frac{\beta_1^{\top}x_1+\rho_{12}\frac{\ln y - \beta_2^{\top}x_2}{\sigma}}{\sqrt{1-\rho_{12}^2}}\right)$ & $\exp\left\{\beta_2^{\top} x_2 + \frac{\sigma^2}{2}\right\}\frac{\Phi\left(\beta_1^{\top}x_1+\sigma\rho_{12}\right)}{\Phi_1}$\\
3 &$\square$&$\blacksquare$&$\square$&$\square$&$\square$&$\square$& $1 - \Phi_2$ & $\frac{1}{\sigma}\phi\left(\frac{y-\beta_2^{\top}x_2}{\sigma}\right)$ & $\beta_2^{\top}x_2+\sigma\frac{\phi_2}{\Phi_2}$\\
4i &$\square$&$\square$&$\blacksquare$&$\square$&$\square$&$\square$& $1 - \Phi_3$ & $\frac{1}{\sigma y}\phi\left(\frac{\ln(\Phi_3 y)-\beta_2^{\top}x_2}{\sigma}\right)\Phi_3$ & $\exp\left\{\beta_2^{\top} x_2 + \frac{\sigma^2}{2}\right\}\frac{1}{\Phi_3}$ \\
4d &$\square$&$\square$&$\blacksquare$&$\square$&$\square$ &$\blacksquare$& $1 - \Phi_3$ & $\frac{1}{\sigma y}\phi\left(\frac{\ln (\Phi_3 y)-\beta_2^{\top}x_2}{\sigma}\right)\Phi\left(\frac{\beta_3^{\top}x_3 + \rho_{23}\frac{\ln (\Phi_3 y) - \beta_2^{\top}x_2}{\sigma}}{\sqrt{1-\rho_{23}^2}}\right)$  & $\exp\left\{\beta_2^{\top} x_2
+ \frac{\sigma^2}{2}\right\}\frac{\Phi\left(\beta_3^{\top}x_3+\sigma\rho_{23}\right)}{\Phi_3^2}$\\
5i &$\blacksquare$&$\blacksquare$&$\square$&$\square$&$\square$&$\square$& $1 - \Phi_1\Phi_2$ & $\frac{1}{\sigma}\phi\left(\frac{y-\beta_2^{\top}x_2}{\sigma}\right)\Phi_1$ & $\beta_2^{\top} x_2 + \sigma\frac{\phi_2}{\Phi_{2}}$\\
5d &$\blacksquare$&$\blacksquare$&$\square$&$\blacksquare$&$\square$&$\square$&$1 - \Phi_{12}$ &
$\frac{1}{\sigma}\phi\left(\frac{y-\beta_2^{\top}x_2}{\sigma}\right)\Phi\left(\frac{\beta_1^{\top}x_1 + \rho_{12}\frac{y-\beta_2^{\top}x_2}{\sigma}}{\sqrt{1-\rho_{12}^2}}\right)$ & $\beta_2^{\top} x_2 + \sigma \frac{\Psi_{2|1}}{\Phi_{12}}$\\
6i &$\blacksquare$&$\square$&$\blacksquare$&$\square$&$\square$&$\square$& $1 - \Phi_1 \Phi_3$ & $\frac{1}{\sigma y}\phi\left(\frac{\ln (\Phi_3 y) -\beta_2^{\top}x_2}{\sigma}\right)\Phi_1\Phi_3$ &$\exp\left\{\beta_2^{\top}x_2+\frac{\sigma^2}{2}\right\}\frac{1}{\Phi_3}$ \\
6d1 &$\blacksquare$&$\square$&$\blacksquare$&$\blacksquare$&$\square$&$\square$& $1 - \Phi_1 \Phi_3$ &  $\frac{1}{\sigma y}\phi\left(\frac{(\ln \Phi_3 y) -\beta_2^{\top}x_2}{\sigma}\right)\Phi\left(\frac{\beta_1^{\top}x_1+\rho_{12}\frac{(\ln y \Phi_3) - \beta_2^{\top}x_2}{\sigma}}{\sqrt{1-\rho_{12}^2}}\right)\Phi_3$&
$\exp\left\{\beta_2^{\top} x_2 + \frac{\sigma^2}{2}\right\}\frac{\Phi\left(\beta_1^{\top}x_1+\sigma\rho_{12}\right)}{\Phi_1\Phi_3}$\\
6d3 &$\blacksquare$&$\square$&$\blacksquare$&$\square$&$\square$ &$\blacksquare$& $1 - \Phi_1 \Phi_3$ & $\frac{1}{\sigma y}\phi\left(\frac{\ln (\Phi_3 y)-\beta_2^{\top}x_2}{\sigma}\right)\Phi_1\Phi\left(\frac{\beta_3^{\top}x_3 + \rho_{23}\frac{\ln (\Phi_3 y) - \beta_2^{\top}x_2}{\sigma}}{\sqrt{1-\rho_{23}^2}}\right)$& $\exp\left\{\beta_2^{\top} x_2
+\frac{\sigma^2}{2}\right\}\frac{\Phi\left(\beta_3^{\top}x_3+\sigma\rho_{23}\right)}{\Phi_3^2}$\\
7i &$\square$&$\blacksquare$&$\blacksquare$&$\square$&$\square$&$\square$& $1 - \Phi_2 \Phi_3$& $\frac{1}{\sigma}\phi\left(\frac{\Phi_3y-\beta_2^{\top}x_2}{\sigma}\right)\Phi_3^2$ & $\frac{\beta_2^{\top} x_2}{\Phi_3}
+\sigma \frac{\phi_2}{\Phi_2 \Phi_3}$\\
7d &$\square$&$\blacksquare$&$\blacksquare$&$\square$&$\square$ &$\blacksquare$& $1 - \Phi_{23}$& $\frac{1}{\sigma}\phi\left(\frac{\Phi_3 y-\beta_2^{\top}x_2}{\sigma}\right)\Phi\left(\frac{\beta_3^{\top}x_3 + \rho_{23}\frac{\Phi_3 y - \beta_2^{\top}x_2}{\sigma}}{\sqrt{1-\rho_{23}^2}}\right)\Phi_3$ & $\frac{\beta_2^{\top} x_2}{\Phi_3}
+ \sigma \frac{\Psi_{2|3}}{\Phi_{23}\Phi_3}$ \\
8i  &$\blacksquare$&$\blacksquare$&$\blacksquare$&$\square$&$\square$&$\square$& $1-\Phi_1\Phi_2\Phi_3$& $\frac{1}{\sigma}\phi\left(\frac{\Phi_3 y-\beta_2^{\top}x_2}{\sigma}\right)\Phi_1\Phi_3^2$& $\frac{\beta_2^{\top} x_2}{\Phi_3} + \sigma \frac{\phi_2}{\Phi_2 \Phi_3}$ \\
8d1 &$\blacksquare$&$\blacksquare$&$\blacksquare$&$\blacksquare$&$\square$&$\square$&$1 - \Phi_{12}\Phi_3$& $\frac{1}{\sigma}\phi\left(\frac{\Phi_3 y-\beta_2^{\top}x_2}{\sigma}\right)\Phi\left(\frac{\beta_1^{\top}x_1 + \rho_{12}\frac{\Phi_3 y - \beta_2^{\top}x_2}{\sigma}}{\sqrt{1-\rho_{12}^2}}\right)\Phi_3^2$ & $\frac{\beta_2^{\top} x_2}{\Phi_3}
+ \sigma \frac{\Psi_{2|1}}{\Phi_{12} \Phi_3}$ \\
8d3 &$\blacksquare$&$\blacksquare$&$\blacksquare$&$\square$&$\square$ &$\blacksquare$&$1 - \Phi_1\Phi_{23}$& $\frac{1}{\sigma}\phi\left(\frac{\Phi_3 y-\beta_2^{\top}x_2}{\sigma}\right)\Phi_1\Phi\left(\frac{\beta_3^{\top}x_3 + \rho_{23}\frac{\Phi_3 y - \beta_2^{\top}x_2}{\sigma}}{\sqrt{1-\rho_{23}^2}}\right)\Phi_3$ & $\frac{\beta_2^{\top} x_2}{\Phi_3}
+ \sigma \frac{\Psi_{2|3}}{\Phi_{23} \Phi_3}$  \\ \hline\hline
\end{tabular}

%% BC <- mhurdle(comics ~ gender + educ + age |  log(incum) +
%%                  I(log(incum)^2) + I(log(incum)^3) + size | 0,
%%                  data = Comics, corr = "d", dist = "bc", method = 'bfgs', print.level = 3)
                 
%% IHS <- mhurdle(comics ~ gender + educ + age |  log(incum) +
%%                  I(log(incum)^2) + I(log(incum)^3) + size | 0,
%%                  data = Comics, corr = "d", dist = "ihs", method = 'bfgs', print.level = 3)                 

                 
%% <<label = onequation>>=
%% Sn <- mhurdle(vacations ~ 0 | linc,# + linc2,# + size + smsa,
%%               Interview, dist = "n", h2 = TRUE, corr = FALSE,
%%               method = "bhhh", print.level = 0)
%% Sbc <- update(Sn, dist = "bc")
%% Sln <- update(Sn, dist = "ln")
%% library(texreg)
%% @ 
%% <<eval = FALSE, echo = TRUE, results = hide>>=
%% texreg(list(Sn, Sln, Sbc), 
%%        custom.model.names = c("normal tobit", "log-normal tobit", "box-cox tobit"), 
%%        caption = "Estimation of one-equation models", 
%%        label = "tab:oneq", pos = "ht")
%% @ 

                
%% <<label = onequation>>=
%% Sn <- mhurdle(tobacco ~ age + educ | linc,# + linc2,# + size + smsa,
%%               Interview, dist = "n", h2 = FALSE, corr = FALSE,
%%               method = "bhhh", print.level = 0)
%% Sbc <- update(Sn, dist = "bc")
%% Sln <- update(Sn, dist = "ln")
%% library(texreg)
%% @ 
%% <<eval = FALSE, echo = TRUE, results = hide>>=
%% texreg(list(Sn, Sln, Sbc), 
%%        custom.model.names = c("normal tobit", "log-normal tobit", "box-cox tobit"), 
%%        caption = "Estimation of one-equation models", 
%%        label = "tab:oneq", pos = "ht", single.row=TRUE)
%% @ 

%% <<label = onequation>>=
%% Sn <- mhurdle(apparel ~ 0 | linc|smsa, Interview, dist = "ln", 
%%               h2 = FALSE, corr = FALSE,      method = "bfgs", print.level = 0)
%% Sbc <- update(Sn, dist = "bc")
%% Sln <- update(Sn, dist = "ln")
%% library(texreg)
%% @ 
%% <<eval = FALSE, echo = TRUE, results = hide>>=
%% texreg(list(Sn, Sln, Sbc), 
%%        custom.model.names = c("normal tobit", "log-normal tobit", "box-cox tobit"), 
%%        caption = "Estimation of one-equation models", 
%%        label = "tab:oneq", pos = "ht", single.row=TRUE)
%% @ 
