Spending and length of stay by tourists flying to the Canary Islands (Spain) using low-cost carriers

Abstract

The AIM of this study is to examine and evaluate differences in expenditure and length of stay between tourists who use low-cost carriers and those who travel with full service providers. We consider the statistical dependence between these variables and propose a bivariate distribution that describes tourist expenditure (continuous variable) and length of stay (discrete variable) in terms of their conditional distributions. Covariates are included to reflect the factors that simultaneously affect both variables. In addition, an empirical analysis is made of data obtained by the Canary Islands Tourist Expenditure Survey. The results obtained show that our model achieves a reasonably good fit and that there are differences between LCC and FSC users regarding both expenditure and length of stay, in the use of nonhotel accommodation, as well as differences in expenditure in the case of repeated visits, and in the length of stay according to the visitors’ age, nationality and travel party size.

Appendix A

A.1 Marginal, conditional distributions and marginal moments

The marginal distributions of X and N are given by

$$ \begin{array}{@{}rcl@{}} X & \sim & \mathcal{G}(\alpha,\gamma),\quad \alpha>0, \beta>\gamma, \gamma > 0, \end{array} $$

(15)

$$ \begin{array}{@{}rcl@{}} N & \sim & \mathcal{S}\mathcal{N}\mathcal{B}\left( \alpha,p=\frac{\gamma}{\beta}\right),\quad \alpha>0, \beta>\gamma, \gamma > 0, \end{array} $$

(16)

where $\mathcal {S}\mathcal {N}{\mathscr{B}}$ refers to the shifted negative binomial distribution. That is,

$$ \begin{array}{@{}rcl@{}} f_{X}(x) &=& \frac{\gamma^{\alpha}}{{{{\varGamma}}}(\alpha)}x^{\alpha-1}\exp(-\gamma x),\quad x>0, \end{array} $$

(17)

$$ \begin{array}{@{}rcl@{}} f_{N}(n) &=& {\alpha+n-2\choose n-1}\left( \frac{\gamma}{\beta}\right)^{\alpha}\left( 1-\frac{\gamma}{\beta}\right)^{n-1}, n=1,2,{\dots} \end{array} $$

(18)

The marginal distribution of N is obtained by integrating (13) with respect to x in the support $(0,+\infty )$ and the marginal distribution of X is calculated by summing (13) with respect to n in the support $\{1,2,\dots \}$.

Furthermore, it is readily apparent that the marginal means and variances are given by

$$ \begin{array}{@{}rcl@{}} E(X) &=& \frac{\alpha}{\gamma},\quad var(X)= \frac{\alpha}{\gamma^{2}},\\ E(N) &=& 1+\frac{\alpha}{\gamma}(\beta-\gamma),\quad var(N) = \frac{\alpha\beta(\beta-\gamma)}{\gamma^{2}}, \end{array} $$

while the cross moment of X and N is

$$ \begin{array}{@{}rcl@{}} E(XN) = \frac{\alpha}{\gamma^{2}}[\beta+\alpha(\beta-\gamma)]. \end{array} $$

Simple calculations provide the covariance, given by

$$ \begin{array}{@{}rcl@{}} cov(X,N) = \frac{\alpha(\beta-\gamma)}{\gamma^{2}}, \end{array} $$

which is always positive.

The conditional distribution of X given N = n is $\mathcal {G}(\sigma (n),\eta (n))$, where

$$ \begin{array}{@{}rcl@{}} \sigma(n) &=& \alpha+n-1,\\ \eta(n) &=& \beta \end{array} $$

and the conditional distribution of N given X = x is $\mathcal {S}\mathcal {P}o(\varphi (x))$, with φ(x) = (β − γ)x. Observe that with the assumption of these two conditional distributions expressions (1) and (2) are guaranteed.

A.2 Estimation of the parameters

Here, we derive estimators based on the moments method and on maximum likelihood for the model with and without covariates, and also provide closed-form expressions for the Fisher information matrix.

A.3 Estimation of the model without covariates

Let us first consider the case of the model with no covariates. If

$$ \begin{array}{@{}rcl@{}} (\tilde x,\tilde n)=\{(x_{1},n_{1}),(x_{2},n_{2}),\dots, (x_{t},n_{t})\} \end{array} $$

is a sample obtained from the distribution (14) and $\bar x=(1/t){\sum }_{i=1}^{t}x_{i}$, $\bar n=(1/t){\sum }_{i=1}^{t}n_{i}$ and $\mu _{12}=(1/t){\sum }_{i=1}^{t}x_{i} n_{i}$ are the corresponding sample moments, some computation provides the estimators based on these sample moments, which are given by

$$ \begin{array}{@{}rcl@{}} \widehat\mu_{1} = \bar x,\quad \widehat\mu_{2}=\bar n,\quad\widehat\gamma=\frac{\widehat\mu_{2}-1}{\widehat \mu_{12}-\widehat \mu_{1} \widehat \mu_{2}}. \end{array} $$

A.4 The score vector and Fisher information matrix

We now consider the maximum likelihood method. Let Θ = (γ,μ₁,μ₂) be the vector of parameters to be estimated. The log-likelihood function is proportional to

$$ \begin{array}{@{}rcl@{}} \ell((\tilde x,\tilde n);{{{\varTheta}}}) & \propto & t\left[\gamma\mu_{1}(\bar x^{\ast}+\log\gamma)-\bar x^{\ast}-\log{{{\varGamma}}}(\gamma\mu_{1})\right]+\sum\limits_{i=1}^{t}n_{i}\log x_{i}\\ &&+t(\bar n-1)\left[\log(\mu_{2}-1)-\log\mu_{1}\right] -\frac{t\bar x(\mu_{2}+\gamma\mu_{1}-1)}{\mu_{1}}, \end{array} $$

where $\bar x^{\ast }=(1/t){\sum }_{i=1}^{t}\log x_{i}$.

Thus, the normal equations which provide the estimators of the parameters are given by

$$ \begin{array}{@{}rcl@{}} \bar x^{\ast}+\log\gamma+1-\psi(\gamma\mu_{1})-\frac{\bar x}{\mu_{1}} &=& 0, \end{array} $$

(19)

$$ \begin{array}{@{}rcl@{}} {\mu_{1}^{2}}\gamma\left[\bar x^{\ast}+\log\gamma-\psi(\gamma\mu_{1})\right]-\mu_{1}(\bar n-1)+\bar x(\mu_{2}-1) &=& 0, \end{array} $$

(20)

$$ \begin{array}{@{}rcl@{}} \mu_{1}(\bar n-1)-\bar x(\mu_{2}-1) &=& 0, \end{array} $$

(21)

where ψ(z) is the digamma function, the logarithmic derivative of the Euler gamma function. Some algebra manipulation provides the maximum likelihood estimators of μ₁ and μ₂, which are given by $\widehat \mu _{1}=\bar x$ and $\widehat \mu _{2}=\bar n$. Finally, the estimator of the parameter γ is the solution of the equation

$$ \begin{array}{@{}rcl@{}} \log\gamma-\psi(\gamma\mu_{1})+\bar x^{\ast}=0, \end{array} $$

which can be solved numerically.

The second partial derivatives are as follows:

$$ \begin{array}{@{}rcl@{}} \frac{\partial^{2} \ell((\tilde x,\tilde n);{{{\varTheta}}})}{\partial\gamma^{2}} &=& \frac{t\mu_{1}}{\gamma}-t{\mu_{1}^{2}}\psi_{1}(\gamma\mu_{1}),\\ \frac{\partial^{2} \ell((\tilde x,\tilde n);{{{\varTheta}}})}{\partial\gamma\partial\mu_{1}} &=& t\left[\bar 1+x^{\ast}+\log\gamma-\psi(\gamma\mu_{1})-\gamma\mu_{1}\psi_{1}(\gamma\mu_{1})\right],\\ \frac{\partial^{2} \ell((\tilde x,\tilde n);{{{\varTheta}}})}{\partial\gamma\mu_{2}} &=& 0,\\ \frac{\partial^{2} \ell((\tilde x,\tilde n);{{{\varTheta}}})}{{\partial\mu_{1}^{2}}} &=& \frac{t(\bar n-1)}{{\mu_{1}^{2}}}-t\gamma^{2}\psi_{1}(\gamma\mu_{1})-\frac{2t\bar x}{{\mu_{1}^{3}}}(\mu_{2}-1),\\ \frac{\partial^{2} \ell((\tilde x,\tilde n);{{{\varTheta}}})}{\partial\mu_{1}\mu_{2}} &=& \frac{t\bar x}{{\mu_{1}^{2}}},\quad \frac{\partial^{2} \ell((\tilde x,\tilde n);{{{\varTheta}}})}{{\partial\mu_{2}^{2}}} = -\frac{t(\bar n-1)}{(\mu_{2}-1)^{2}}, \end{array} $$

where ψ₁(⋅) is the first derivative of the digamma function.

The elements of the Fisher information matrix, $\mathcal {J}(\widehat {{{\varTheta }}})$, are therefore

$$ \begin{array}{@{}rcl@{}} \mathcal{J}_{11}(\widehat{{{\varTheta}}}) &=& \left.E\left( -\frac{\partial^{2} \ell((\tilde x,\tilde n);{{{\varTheta}}})}{\partial\gamma^{2}}\right)\right|_{{{{\varTheta}}}=\widehat{{{\varTheta}}}}= -\frac{t\widehat\mu_{1}}{\widehat\gamma}+t{\widehat\mu_{1}^{2}}\psi_{1}(\widehat\gamma\widehat\mu_{1}),\\ \mathcal{J}_{12}(\widehat{{{\varTheta}}}) & = & \left.E\left( -\frac{\partial^{2} \ell((\tilde x,\tilde n);{{{\varTheta}}})}{\partial\widehat\gamma\partial\widehat\mu_{1}}\right)\right|_{{{{\varTheta}}}=\widehat{{{\varTheta}}}}=-\sum\limits_{i=1}^{t}\log x_{i}-t\left[(1+\log\widehat\gamma-\psi(\widehat\gamma\widehat\mu_{1})\right.\\ &&\left.-\widehat\gamma\widehat\mu_{1}\psi_{1}(\widehat\gamma\widehat\mu_{1})\right],\\ \mathcal{J}_{13}(\widehat{{{\varTheta}}}) &=& \left.E\left( -\frac{\partial^{2} \ell((\tilde x,\tilde n);{{{\varTheta}}})}{\partial\widehat\gamma\widehat\mu_{2}}\right)\right|_{{{{\varTheta}}}=\widehat{{{\varTheta}}}}= 0,\\ \mathcal{J}_{22}(\widehat{{{\varTheta}}}) &=& \left.E\left( -\frac{\partial^{2} \ell((\tilde x,\tilde n);{{{\varTheta}}})}{{\partial\widehat\mu_{1}^{2}}}\right)\right|_{{{{\varTheta}}}=\widehat{{{\varTheta}}}}= -\frac{t(\widehat\mu_{2}-1)}{{\widehat\mu_{1}^{2}}}+t\widehat\gamma^{2}\psi_{1}(\widehat\gamma\widehat\mu_{1})\\ &&+\frac{2t}{{\widehat\mu_{1}^{2}}}(\widehat\mu_{2}-1),\\ \mathcal{J}_{23}(\widehat{{{\varTheta}}}) &=& \left.E\left( -\frac{\partial^{2} \ell((\tilde x,\tilde n);{{{\varTheta}}})}{\partial\widehat\mu_{1}\widehat\mu_{2}}\right)\right|_{{{{\varTheta}}}=\widehat{{{\varTheta}}}}= -\frac{t}{\widehat\mu_{1}},\\ \mathcal{J}_{33}(\widehat{{{\varTheta}}}) &=& \left.E\left( -\frac{\partial^{2} \ell((\tilde x,\tilde n);{{{\varTheta}}})}{{\partial\widehat\mu_{2}^{2}}}\right)\right|_{{{{\varTheta}}}=\widehat{{{\varTheta}}}}=\frac{t}{(\widehat\mu_{2}-1)}. \end{array} $$

Here, $\widehat {{{\varTheta }}}$ represents the maximum likelihood of Θ. Observe that the analytic expression for $E({\sum }_{i=1}^{t}\log x_{i})$ is not feasible. For large t, for computational purposes, this is evaluated by ignoring the expectation operator and replacing it by ${\sum }_{i=1}^{t} \log x_{i}$. The asymptotic variance-covariance matrix of $\widehat {{{\varTheta }}}$ is obtained by inverting the information matrix.

A.5 Estimation of the model with covariates

For the sake of simplicity, we assume η = δ and write μ_1i = μ_1i(δ) and μ_2i = μ_2i(δ). Let Θ = (γ,δ). The log-likelihood is then proportional to

$$ \begin{array}{@{}rcl@{}} \ell((\tilde x,\tilde n);{{{\varTheta}}}) &\propto& \gamma\sum\limits_{i=1}^{t} \mu_{1i}\log(\gamma x_{i})-\sum\limits_{i=1}^{t}\log{{{\varGamma}}}(\gamma\mu_{1i})\\ &&+ \sum\limits_{i=1}^{t}(n_{i}-1)\left[\log x_{i}+\log(\mu_{2i}-1)-\log\mu_{1i}\right]\\ &&-\sum\limits_{i=1}^{t}\frac{x_{i}}{\mu_{1i}}(\mu_{2i}+\gamma\mu_{1i}-1). \end{array} $$

Thus, the normal equations, for $i=1,\dots ,t$, are given by

$$ \begin{array}{@{}rcl@{}} \frac{\partial \ell((\tilde x,\tilde n);{{{\varTheta}}})}{\partial \gamma} &=& \sum\limits_{i=1}^{t} \mu_{1i}\left[1+\log(\gamma x_{i})-\psi(\gamma \mu_{1i})\right]-t\bar x=0,\\ \frac{\partial \ell((\tilde x,\tilde n);{{{\varTheta}}})}{\partial \delta_{j}} &=& \gamma\sum\limits_{i=1}^{t} y_{ij}\mu_{1i}\left[\log(\gamma x_{i})-\psi(\gamma\mu_{1i})\right]\\ &&+\sum\limits_{i=1}^{t}(n_{i}-1)\left[\frac{1}{\mu_{2i}-1}\frac{\partial\mu_{2i}}{\partial \delta_{j}}-\frac{1}{\mu_{1i}}\frac{\partial\mu_{1i}}{\partial\delta_{j}}\right]\\ &&-\sum\limits_{i=1}^{t}\frac{x_{i}}{\mu_{1i}^{2}}\left[\mu_{1i}\frac{\partial\mu_{2i}}{\partial \delta_{j}}+(1-\mu_{2i})\frac{\partial\mu_{1i}}{\partial\delta_{j}}\right]=0, \end{array} $$

where $j=1,\dots ,k$, $\frac {\partial \mu _{1i}}{\partial \delta _{j}} = y_{ij}\mu _{1i}$ and $\frac {\partial \mu _{2i}}{\partial \delta _{j}} = z_{ij}(\mu _{2i}-1)$. Finally, after computing the second partial derivatives we obtain the elements of the Fisher information matrix, as follows:

$$ \begin{array}{@{}rcl@{}} \mathcal{J}_{11}(\widehat{{{\varTheta}}}) &=& \left.E\left( -\frac{\partial^{2} \ell((\tilde x,\tilde n);{{{\varTheta}}})}{\partial\gamma^{2}}\right)\right|_{{{{\varTheta}}}=\widehat{{{\varTheta}}}}=-\sum\limits_{i=1}^{t} \widehat\mu_{1i}\left[\frac{1}{\widehat\gamma}-\widehat\mu_{1i}\psi_{1}(\widehat\gamma\widehat\mu_{1i})\right],\\ \mathcal{J}_{12}(\widehat{{{\varTheta}}}) &=& \left.E\left( -\frac{\partial^{2} \ell((\tilde x,\tilde n);{{{\varTheta}}})}{\partial\gamma\partial\delta_{j}}\right)\right|_{{{{\varTheta}}}=\widehat{{{\varTheta}}}}=-\sum\limits_{i=1}^{t} y_{ij}\widehat\mu_{1i}\left[1+\log(\widehat\gamma x_{i})\right.\\ &&\left.-\psi(\widehat\gamma\widehat\mu_{1i})-\widehat\gamma\widehat\mu_{1i} \psi_{1}(\widehat\gamma\widehat\mu_{1i})\right],\\ \mathcal{J}_{22}(\widehat{{{\varTheta}}}) &=& \left.E\left( -\frac{\partial^{2} \ell((\tilde x,\tilde n);{{{\varTheta}}})}{{\partial\delta_{j}^{2}}}\right)\right|_{{{{\varTheta}}}=\widehat{{{\varTheta}}}}=-\sum\limits_{i=1}^{t} y_{ij}^{2}\mu_{1i}\widehat\gamma\left[\log(\widehat\gamma x_{i})\right.\\ &&\left.-\psi(\widehat\gamma \widehat\mu_{1i})\right]+y_{ij}^{2}\widehat\mu_{1i}^{2}\left[-\widehat\gamma^{2}\psi_{1}(\widehat\gamma \widehat\mu_{1i})+\frac{1-\mu_{2i}}{\mu_{1i}}\right]+y_{ij}^{2}(\mu_{2i}-1),\\ \mathcal{J}_{23}(\widehat{{{\varTheta}}}) &=& \left.E\left( -\frac{\partial^{2} \ell((\tilde x,\tilde n);{{{\varTheta}}})}{\partial\delta_{j}\partial\delta_{l}}\right)\right|_{{{{\varTheta}}}=\widehat{{{\varTheta}}}}=-\sum\limits_{i=1}^{t} y_{ij}z_{il}\widehat\gamma\widehat\mu_{1i}\left[\log(\widehat\gamma x_{i})-\psi(\widehat\gamma\widehat\mu_{1i})\right]\\ &&+y_{ij}z_{il}\widehat\mu_{1i}^{2}\left[-\widehat\gamma^{2}\psi_{1}(\widehat\gamma\widehat\mu_{1i})+ \frac{1-\widehat\mu_{2i}}{\widehat\mu_{1i}^{2}}\right]+z_{il}^{2}(\widehat\mu_{2i}-1), j\neq l. \end{array} $$

Standard errors can be obtained conventionally from the inverse of the matrix.

A.6 Approximation of the digamma function

A practical approach to the digamma function is given by the following expression, which is well known in the statistical literature,

$$ \begin{array}{@{}rcl@{}} \log\left( {{{\varGamma}}}(z)\right)\approx \frac{1}{2}\log(2\pi)+\left( z-\frac{1}{2}\right) \log(z)-z+\frac{z}{2}\log\left( z\sinh\left( \frac{1}{z}\right)\right). \end{array} $$

(22)