Normal distribution

The density of normal distribution 𝒩(μ_t, σ) is: where σ is the standard deviation of the error term. This PDF has a very well-known bell shape:

alm() with Normal distribution (distribution="dnorm") is equivalent to lm() function from stats package and returns roughly the same estimates of parameters, so if you are concerned with the time of calculation, I would recommend reverting to lm().

Maximising the likelihood of the model is equivalent to the estimation of the basic linear regression using Least Squares method: where ϵ_t ∼ 𝒩(0, σ²).

The variance σ² is estimated in alm() based on likelihood: where T is the sample size. Its square root (standard deviation) is used in the calculations of dnorm() function, and the value is then return via scale variable. This value does not have bias correction. However the sigma() method applied to the resulting model, returns the bias corrected version of standard deviation. And vcov(), confint(), summary() and predict() rely on the value extracted by sigma().

μ_t is returned as is in mu variable, and the fitted values are set equivalent to mu.

In order to produce confidence intervals for the mean (predict(model, newdata, interval="confidence")) the conditional variance of the model is calculated using: where V(B) is the covariance matrix of the parameters returned by the function vcov. This variance is then used for the construction of the confidence intervals of a necessary level α using the distribution of Student: where $\tau_{df,\frac{1+\alpha}{2}}$ is the upper ${\frac{1+\alpha}{2}}$-th quantile of the Student’s distribution with df degrees of freedom (e.g. with α = 0.95 it will be 0.975-th quantile, which, for example, for 100 degrees of freedom will be  ≈ 1.984).

Similarly for the prediction intervals (predict(model, newdata, interval="prediction")) the conditional variance of the y_t is calculated: where s² is the bias-corrected variance of the error term, calculated using: where k is the number of estimated parameters (including the variance itself). This value is then used for the construction of the prediction intervals of a specify level, also using the distribution of Student, in a similar manner as with the confidence intervals.

Laplace distribution

Laplace distribution has some similarities with the Normal one: where s is the scale parameter, which, when estimated using likelihood, is equal to the mean absolute error: So maximising the likelihood is equivalent to estimating the linear regression via the minimisation of s . So when estimating a model via minimising s, the assumption imposed on the error term is ϵ_t ∼ ℒ𝒶𝓅𝓁𝒶𝒸ℯ(0, s). The main difference of Laplace from Normal distribution is its fatter tails, the PDF has the following shape:

alm() function with distribution="dlaplace" returns mu equal to μ_t and the fitted values equal to mu. s is returned in the scale variable. The prediction intervals are derived from the quantiles of Laplace distribution after transforming the conditional variance into the conditional scale parameter s using the connection between the two in Laplace distribution: where σ² is substituted either by the conditional variance of μ_t or y_t.

The kurtosis of Laplace distribution is 6, making it suitable for modelling rarely occurring events.

Asymmetric Laplace distribution

Asymmetric Laplace distribution can be considered as a two Laplace distributions with different parameters s for left and right side. There are several ways to summarise the probability density function, the one used in alm() relies on the asymmetry parameter α (Yu and Zhang 2005): where s is the scale parameter, α is the skewness parameter and I(y_t ≤ μ_t) is the indicator function, which is equal to one, when the condition is satisfied and to zero otherwise. The scale parameter s estimated using likelihood is equal to the quantile loss: Thus maximising the likelihood is equivalent to estimating the linear regression via the minimisation of α quantile, making this equivalent to quantile regression. So quantile regression models assume indirectly that the error term is ϵ_t ∼ 𝒜ℒ𝒶𝓅𝓁𝒶𝒸ℯ(0, s, α) (Geraci and Bottai 2007). The advantage of using alm() in this case is in having the full distribution, which allows to do all the fancy things you can do when you have likelihood.

Graphically, the PDF of asymmetric Laplace is:

In case of α = 0.5 the function reverts to the symmetric Laplace where $s=\frac{1}{2}\text{MAE}$.

alm() function with distribution="dalaplace" accepts an additional parameter alpha in ellipsis, which defines the quantile α. If it is not provided, then the function will estimated it maximising the likelihood and return it as the first coefficient. alm() returns mu equal to μ_t and the fitted values equal to mu. s is returned in the scale variable. The parameter α is returned in the variable other of the final model. The prediction intervals are produced using qalaplace() function. In order to find the values of s for the holdout the following connection between the variance of the variable and the scale in Asymmetric Laplace distribution is used: where σ² is substituted either by the conditional variance of μ_t or y_t.

NOTE: in order for the Asymmetric Laplace to work well, you might need to have large samples. This is inherited from the pinball score of the quantile regression. If you fit the model on 40 observations with α = 0.05, you will only have 2 observations below the line, which does not help very much with the fit. Similarly, the covariance matrix, produced via the Hessian might not be adequate in this situation (because there is not enough variability in the data due to extreme value of α). The latter can be partially addressed by using bootstrap, but do not expect miracles on small samples.

S distribution

The S distribution has the following density function: where s is the scale parameter. If estimated via maximum likelihood, the scale parameter is equal to: which corresponds to the minimisation of a half of “Mean Root Absolute Error” or “Half Absolute Moment”.

S distribution has a kurtosis of 25.2, which makes it a “severe excess” distribution (thus the name). It might be useful in cases of randomly occurring incidents and extreme values (Black Swans?). Here how the PDF looks:

alm() function with distribution="ds" returns μ_t in the same variables mu and fitted.values, and s in the scale variable. Similarly to the previous functions, the prediction intervals are based on the qs() function from greybox package and use the connection between the scale and the variance: where once again σ² is substituted either by the conditional variance of μ_t or y_t.

Generalised Normal distribution

The Generalised Normal distribution is a generalisation, which has Normal, Laplace and S as special cases. It has the following density function: where s is the scale and β is the shape parameters. If estimated via maximum likelihood, the scale parameter is equal to: In the special cases, this becomes either $\sqrt{2}\times$RMSE (β = 2), or MAE (β = 1) or a half of HAM (β = 0.5). It is important to note that although in case of β = 2, the distribution becomes equivalent to Normal, the scale of it will differ from the σ (this follows directly from the formula above). The relations between the two is: s² = 2σ².

The kurtosis of Generalised Normal distribution is determined by β and is equal to $\frac{\Gamma(5/\beta)\Gamma(1/\beta)}{\Gamma(3/\beta)^2}$.

alm() function with distribution="dgnorm" returns μ_t in the same variables mu and fitted.values, s in the scale variable and β in other$beta. Note that if beta is not provided in the function, then it will estimate it. However, the estimates of β are known not to be consistent and asymptotically normal if it is less than 2. So, use with care! As for the intervals, they are based on the qgnorm() function from greybox package and use the connection between the scale and the variance: where once again σ² is substituted either by the conditional variance of μ_t or y_t, depending on what type of interval is needed.

Logistic distribution

The density function of Logistic distribution is: where s is the scale parameter, which is estimated in alm() based on the connection between the parameter and the variance in the logistic distribution: Once again the maximisation of implies the estimation of the linear model , where ϵ_t ∼ ℒℴℊ𝒾𝓈𝓉𝒾𝒸(0, s).

Logistic is considered a fat tailed distribution, but its tails are not as fat as in Laplace. Kurtosis of standard Logistic is 4.2.

alm() function with distribution="dlogis" returns μ_t in mu and in fitted.values variables, and s in the scale variable. Similar to Laplace distribution, the prediction intervals use the connection between the variance and scale, and rely on the qlogis function.

Student t distribution

The Student t distribution has a difficult density function: where ν is the number of degrees of freedom, which can also be considered as the scale parameter of the distribution. It has the following connection with the in-sample variance of the error (but only for the case, when ν > 2):

Kurtosis of Student t distribution depends on the value of ν, and for the cases of ν > 4 is equal to $\frac{6}{\nu-4}$. When the μ → ∞, the distribution converges to the normal.

alm() function with distribution="dt" estimates the parameters of the model along with the ν (if it is not provided by the user as a nu parameter) and returns μ_t in the variables mu and fitted.values, and ν in the scale variable. Both prediction and confidence intervals use qt() function from stats package and rely on the conventional number of degrees of freedom T − k. The intervals are constructed similarly to how it is done in Normal distribution (based on qt() function).

An example of application

In order to see how this works, we will create the following data:

set.seed(41, kind="L'Ecuyer-CMRG")
xreg <- cbind(rnorm(200,10,3),rnorm(200,50,5))
xreg <- cbind(500+0.5*xreg[,1]-0.75*xreg[,2]+rs(200,0,3),xreg,rnorm(200,300,10))
colnames(xreg) <- c("y","x1","x2","Noise")

inSample <- xreg[1:180,]
outSample <- xreg[-c(1:180),]

ALM can be run either with data frame or with matrix. Here’s an example with normal distribution and several levels for the construction of prediction interval:

ourModel <- alm(y~x1+x2, data=inSample, distribution="dnorm")
summary(ourModel)
#> Response variable: y
#> Distribution used in the estimation: Normal
#> Loss function used in estimation: likelihood
#> Coefficients:
#>             Estimate Std. Error Lower 2.5% Upper 97.5%  
#> (Intercept) 383.9826    69.5496   246.7240    521.2412 *
#> x1            0.0055     2.2403    -4.4159      4.4269  
#> x2            1.6701     1.2779    -0.8519      4.1920  
#> 
#> Error standard deviation: 88.1329
#> Sample size: 180
#> Number of estimated parameters: 4
#> Number of degrees of freedom: 176
#> Information criteria:
#>      AIC     AICc      BIC     BICc 
#> 2127.157 2127.386 2139.929 2140.523
plot(predict(ourModel,outSample,interval="p",level=c(0.9,0.95)))

And here’s an example with Asymmetric Laplace and predefined α = 0.95:

ourModel <- alm(y~x1+x2, data=inSample, distribution="dalaplace",alpha=0.95)
summary(ourModel)
#> Warning: Choleski decomposition of hessian failed, so we had to revert to the simple inversion.
#> The estimate of the covariance matrix of parameters might be inaccurate.
#> Warning: Sorry, but the hessian is singular, so we could not invert it.
#> Switching to bootstrap of covariance matrix of parameters.
#> Warning in log(yTransformed): NaNs produced
#> Warning in log(ySorted): NaNs produced
#> Response variable: y
#> Distribution used in the estimation: Asymmetric Laplace with alpha=0.95
#> Loss function used in estimation: likelihood
#> Coefficients:
#>             Estimate Std. Error Lower 2.5% Upper 97.5%  
#> (Intercept) 566.7139    41.3321   485.1437    648.2842 *
#> x1            0.0051     0.0038    -0.0023      0.0125  
#> x2            0.8491     0.7832    -0.6965      2.3947  
#> 
#> Error standard deviation: 167.8399
#> Sample size: 180
#> Number of estimated parameters: 4
#> Number of degrees of freedom: 176
#> Information criteria:
#>      AIC     AICc      BIC     BICc 
#> 2348.783 2349.012 2361.555 2362.148
plot(predict(ourModel,outSample))

Box-Cox Normal distribution

Box-Cox Normal distribution used in the greybox package is defined as a distribution that becomes normal after the Box-Cox transformation. This means that if $x=\frac{y^\lambda+1}{\lambda}$ and x ∼ 𝒩(μ, σ²), then y ∼ BC𝒩(μ, σ²). The density function of the Box-Cox Normal distribution in this case is: where the variance estimated using likelihood is: Depending on the value of λ, we will get different shapes of the density function:

When λ = 0 the distribution transforms to the Log-Normal one.

Estimating the model with Box-Cox Normal distribution is equivalent to estimating the parameters of a linear model after the Box-Cox transform: where ϵ_t ∼ 𝒩(0, σ²) or:

alm() with distribution="dbcnorm" does not transform the provided data and estimates the density directly using dbcnorm() function from greybox with the estimated mean μ_t and the variance . The μ_t is returned in the variable mu, the σ² is in the variable scale, while the fitted.values contains the exponent of μ_t, which, given the connection between the Normal and Box-Cox Normal distributions, corresponds to median of distribution rather than mean. Finally, resid() method returns $e_t = \frac{y_t^{\lambda}-1}{\lambda} - \mu_t$. The lambda parameter can be provided by the user via the lambdaBC in ellipsis.

Folded Normal distribution

Folded Normal distribution is obtained when the absolute value of normally distributed variable is taken: if x ∼ 𝒩(μ, σ²), then |x| ∼ Fold𝒩(μ, σ²). The density function is: which can be graphically represented as:

Conditional mean and variance of Folded Normal are estimated in alm() (with distribution="dfnorm") similarly to how this is done for Normal distribution. They are returned in the variables mu and scale respectively. In order to produce the fitted value (which is returned in fitted.values), the following correction is done: where Φ(⋅) is the CDF of Normal distribution.

The model that is assumed in the case of Folded Normal distribution can be summarised as:

The conditional variance of the forecasts is calculated based on the elements of vcov() (as in all the other functions), the predicted values are corrected in the same way as the fitted values , and the prediction intervals are generated from the qfnorm() function of greybox package. As for the residuals, resid() method returns e_t = y_t − μ_t.

Rectified Normal distribution

Rectified Normal distribution is obtained when all the negative values of normally distributed variable are set to zero: if x ∼ 𝒩(μ, σ), then y = max (0, x) ∼ Rect𝒩(μ, σ). The density function is:

where Φ_x(0, μ, σ) is the CDF and ϕ_x(y_t, μ, σ) is the PDF of the Normal distribution. This can be graphically represented as:

This distribution can be useful in modelling intermittent demand, when the demand sizes are not integer.

Conditional location and scale of Rectified Normal are estimated in alm() (with distribution="drectnorm") similarly to how this is done for Normal distribution. They are returned in the variables mu and scale respectively. In order to produce the fitted value (which is returned in fitted.values), the following formula is used:

The model that is assumed in the case of Rectified Normal distribution is:

The conditional variance of the forecasts is calculated based on the elements of vcov() (as in all the other functions), the predicted values are corrected in the same way as the fitted values , and the prediction intervals are generated from the qrectnorm() function of greybox package. As for the residuals, resid() method returns e_t = y_t − μ_t.

Log-Normal distribution

Log-Normal distribution appears when a normally distributed variable is exponentiated. This means that if x ∼ 𝒩(μ, σ²), then exp x ∼ log𝒩(μ, σ²). The density function of Log-Normal distribution is: where the variance estimated using likelihood is: The PDF has the following shape:

Estimating the model with Log-Normal distribution is equivalent to estimating the parameters of log-linear model: where ϵ_t ∼ 𝒩(0, σ²) or:

alm() with distribution="dlnorm" does not transform the provided data and estimates the density directly using dlnorm() function with the estimated mean μ_t and the variance . If you need a log-log model, then you would need to take logarithms of the external variables. The μ_t is returned in the variable mu, the σ² is in the variable scale, while the fitted.values contains the exponent of μ_t, which, given the connection between the Normal and Log-Normal distributions, corresponds to median of distribution rather than mean. Finally, resid() method returns e_t = log y_t − μ_t.

Inverse Gaussian distribution

Inverse Gaussian distribution is an interesting distribution, which is defined for positive values only and has some properties similar to the properties of the Normal distribution. It has two parameters: location μ_t and scale ϕ (aka “dispersion”). There are different ways to parameterise this distribution, we use the dispersion-based one. The important thing that distinguishes the implementation in alm() from the one in glm() or in any other function is that we assume that the model has the following form: and that ϵ_t ∼ ℐ𝒢(1, ϕ). This means that $y_t \sim \mathcal{IG}\left(\mu_t, \frac{\phi}{\mu_t} \right)$, implying that the dispersion of the model changes together with the conditional expectation. The density function for the error term in this case is: where the dispersion parameter is estimated via maximising the likelihood and is calculated using: Note that in our formulation μ_t = exp (x_t′B), so that the mean is always positive. This implies that we deal with a pure multiplicative model. In addition, we assume that μ_t is just a scale for the distribution, otherwise y_t would not follow the Inverse Gaussian distribution. The density function has following shapes depending on the values of parameters:

alm() with distribution="dinvgauss" estimates the density for y_t using dinvgauss() function from statmod package. The μ_t is returned in the variables mu and fitted.values, the dispersion ϕ is in the variable scale. resid() method returns $e_t = \frac{y_t}{\mu_t}$. Finally, the prediction interval for the regression model are generated using qinvgauss() function from the statmod package.

Gamma distribution

Another popular distribution, defined for positive values only is called “Gamma”. It is parametrised via the shape k and scale σ² and has closed forms for mean and variance: E(x) = kσ², V(x) = kσ⁴.

The important thing that distinguishes the implementation in alm() from the one in glm() or in any other function is that we assume that the model has the following form (similar to the Inverse Gaussian model in alm): and that ϵ_t ∼ Γ(σ⁻², σ²), implying that E(ϵ_t) = kσ² = 1 and V(ϵ_t) = σ². This means that y_t ∼ Γ(σ⁻², σ²μ_t), meaning that the variance of the model changes together with the conditional expectation. The density function for the error term in this case is: where the scale parameter σ² can be estimated via the method of moments based on its relation to the variance: Note that in our formulation μ_t = exp (x_t′B), so that the mean is always positive, which implies that we deal with a pure multiplicative model. In addition, we assume that μ_t is just a scale for the distribution, otherwise y_t would not follow Gamma distribution. All of this makes the model restrictive, but arguably reasonable - otherwise the mean of the distribution might behave uncontrollably.

The density function has following shapes depending on the values of parameters:

alm() with distribution="dgamma" estimates the density for y_t using dgamma() function from stats package. The μ_t is returned in the variables mu and fitted.values, the scale σ² is in the variable scale. resid() method returns $e_t = \frac{y_t}{\mu_t}$. Finally, the prediction interval for the regression model are generated using qgamma() function from the stats package.

Exponential distribution

One peculiar and very specific distribution, which can also be used in modelling is Exponential distribution. It only has one parameter, λ, which regulates both mean and variance: It might be useful in cases, when one wants to model inter-arrival times.

The implementation in alm() relies on the model, similar to the Inverse Gaussian and Gamma models: where ϵ_t ∼ Exp(1), implying that E(ϵ_t) = V(ϵ_t) = 1. This is a very restrictive model, which only works in some special cases. If for some reason the variance and mean are not equal to one in the empirical distribution, then the Exponential one would not be appropriate. But in general the model formulated as above implies that $y_t \sim \mathrm{Exp}\left( \frac{1}{\mu_t} \right)$, meaning that the variance of the model changes together with the conditional expectation. The density function for the error term in this case is: Note that in our formulation μ_t = exp (x_t′B), so that the mean is always positive, which implies that we deal with a pure multiplicative model. In addition, we assume that μ_t is just a scale for the distribution, otherwise y_t would not follow Exponential distribution.

The density function has the following shapes depending on the values of the expectation:

alm() with distribution="dexp" estimates the density for y_t using dexp() function from stats package. The μ_t is returned in the variables mu and fitted.values, the scale is assumed to be equal to one. resid() method returns $e_t = \frac{y_t}{\mu_t}$. Finally, the prediction interval for the regression model are generated using qexp() function from the stats package.

NOTE that if the assumption of E(ϵ_t) = V(ϵ_t) = 1 does not hold, the model will produce unreasonable quantiles.

Log-Laplace distribution

This is based on the exponent of Laplace distribution, which means that the PDF in this case is: The model implemented in the package has similarity with Log-Normal distribution. The MLE scale is: The density function of Log-Laplace has the following shapes:

Estimating the model with Log-Laplace distribution is equivalent to estimating the parameters of log-linear model: where ϵ_t ∼ ℒ𝒶𝓅𝓁𝒶𝒸ℯ(0, σ²). This distribution might be useful if the data has a strong skewness (larger than in case of Log-Normal distribution).

alm() with distribution="dllaplace" uses dlaplace() function with the logarithm of actual values, estimated mean μ_t and the scale . The μ_t is returned in the variable mu, the s is in the variable scale, while the fitted.values contains the exponent of μ_t, which corresponds to median of distribution rather than mean. Finally, resid() method returns e_t = log y_t − μ_t.

Log-S distribution

This is based on the exponent of S distribution, giving the PDF: The model implemented in the package has similarity with Log-Normal and Log-Laplace distributions. The MLE scale is: The shape of the density function of Log-S is similar to Log-Laplace but with even more extreme values:

Estimating the model with Log-S distribution is equivalent to estimating the parameters of log-linear model: where ϵ_t ∼ 𝒮(0, σ²). This distribution can be used for sever seldom right tail cases.

alm() with distribution="dls" uses ds() function with the logarithm of actual values, estimated mean μ_t and the scale . The μ_t is returned in the variable mu, the s is in the variable scale, while the fitted.values contains the exponent of μ_t, which corresponds to median of distribution rather than mean. Finally, resid() method returns e_t = log y_t − μ_t.

Log-Generalised Normal distribution

This is based on the exponent of Generalised Normal distribution, giving the PDF: The model implemented in the package has similarity with Log-Normal, Log-Laplace and Log-S distributions. The MLE scale is: The shapes of the distribution depend on the value of parameters, giving it in some cases very long right tail:

Estimating the model with Log-Generalised Normal distribution is equivalent to estimating the parameters of log-linear model: where ϵ_t ∼ 𝒢𝒩(0, s, β).

alm() with distribution="dlgnorm" uses the dgnorm() function from greybox package with the logarithm of actual values, estimated mean μ_t, the scale and either provided or estimated shape parameter β. The μ_t is returned in the variable mu, the s is in the variable scale and β is in other$beta, while the fitted.values contains the exponent of μ_t, which corresponds to median of distribution rather than mean. Finally, resid() method returns e_t = log y_t − μ_t.

Logit-normal distribution

A random variable follows Logit-normal distribution if its logistic transform follows normal distribution: where y ∈ (0, 1), y ∼ logit𝒩(μ, σ²) and z ∼ 𝒩(μ, σ²). The bounds are not supported, because the variable z becomes infinite. The density function of y is: which has the following shapes: Depending on the values of location and scale, the distribution can be either unimodal or bimodal and can be positively or negatively skewed. Because of its connection with normal distribution, the logit-normal has formulae for density, cumulative and quantile functions. However, the moment generation function does not have a closed form.

The scale of the distribution can be estimated via the maximisation of likelihood and has some similarities with the scale in Log-Normal distribution:

Estimating the model with Log-Normal distribution is equivalent to estimating the parameters of logit-linear model: where ϵ_t ∼ 𝒩(0, σ²) or: where $\mathrm{logit}^{-1}(z)=\frac{\exp(z)}{1+\exp(z)}$ is the inverse logistic transform.

alm() with distribution="dlogitnorm" does not transform the provided data and estimates the density directly using dlogitnorm() function from greybox package with the estimated mean μ_t and the variance . The μ_t is returned in the variable mu, the σ² is in the variable scale, while the fitted.values contains the inverse logistic transform of μ_t, which, given the connection between the Normal and Logit-Normal distributions, corresponds to median of distribution rather than mean. Finally, resid() method returns e_t = logit(y_t) − μ_t.

Beta distribution

Beta distribution is a distribution for a continuous variable that is defined on the interval of (0, 1). Note that the bounds are not included here, because the probability density function is not well defined on them. If the provided data contains either zeroes or ones, the function will modify the values using: and it will warn the user about this modification. This correction makes sure that there are no boundary values in the data, and it is quite artificial and needed for estimation purposes only.

The density function of Beta distribution has the form: where α_t is the first shape parameter and β_t is the second one. Note indices for the both shape parameters. This is what makes the alm() implementation of Beta distribution different from any other. We assume that both of them have underlying deterministic models, so that: and where A and B are the vectors of parameters for the respective shape variables. This allows the function to model any shapes depending on the values of exogenous variables. The conditional expectation of the model is calculated using: while the conditional variance is: Beta distribution has shapes similar to the ones of Logit-Normal one, but with shape parameters regulating respectively the left and right tails of the distribution:

alm() function with distribution="dbeta" returns ŷ_t in the variables mu and fitted.values, and V(y_t) in the scale variable. The shape parameters are returned in the respective variables other$shape1 and other$shape2. You will notice that the output of the model contains twice more parameters than the number of variables in the model. This is because of the estimation of two models: α_t and β_t - instead of one.

Respectively, when predict() function is used for the alm model with Beta distribution, the two models are used in order to produce predicted values for α_t and β_t. After that the conditional mean mu and conditional variance variances are produced using the formulae above. The prediction intervals are generated using qbeta function with the provided shape parameters for the holdout. As for the confidence intervals, they are produced assuming normality for the parameters of the model and using the estimate of the variance of the mean based on the variances (which is weird and probably wrong).

Poisson distribution

Poisson distribution used in ALM has the following standard probability mass function (PMF): where λ_t = μ_t = σ_t² = exp (x_t′B). As it can be noticed, here we assume that the variance of the model varies in time and depends on the values of the exogenous variables, which is a specific case of heteroscedasticity. The exponent of x_t′B is needed in order to avoid the negative values in λ_t.

Here are several examples of the PMF of Poisson with different values of parameters λ:

alm() with distribution="dpois" returns mu, fitted.values and scale equal to λ_t. The quantiles of distribution in predict() method are generated using qpois() function from stats package. Finally, the returned residuals correspond to y_t − μ_t, which is not really helpful or meaningful.

Negative Binomial distribution

Negative Binomial distribution implemented in alm() is parameterised in terms of mean and variance: where μ_t = exp (x_t′B) and σ² is estimated separately in the optimisation process. These values are then used in the dnbinom() function in order to calculate the log-likelihood based on the distribution function.

Here are some examples of PMF of Negative Binomial distribution with different sizes and probabilities:

alm() with distribution="dnbinom" returns μ_t in mu and fitted.values and σ² in scale. The prediction intervals are produces using qnbinom() function. Similarly to Poisson distribution, resid() method returns y_t − μ_t. The user can also provide size parameter in ellipsis if it is reasonable to assume that it is known.

Binomial distribution

The PMF of Binomial distribution is written as: where n is the size parameter, μ_t = exp (x_t′B) $p_t = \frac{1}{1+\mu_t}$ and E(y_t) = n × p_t. The size parameter is always known and can be calculated as a number of unique values in y minus one. So, if the data takes one of the three values: 0, 1 and 2, the size will be n = 2. The values of p_t and n are then used in the dbinom() function to calculate the log-likelihood.

Visually, the PMF of the Binomial distribution has the following shapes with different values of probability and size:

alm() with distribution="dbinom" returns μ_t in mu and E(y_t) in fitted.values, the scale and other$size both contain the same value of n. The prediction intervals are produces using qbinom() function. resid() method returns y_t − E(y_t).

Geometric distribution

There is also Geometric distribution implemented in alm(). It has the following PMF: where $p_t = \frac{1}{1+\exp(x_t' B)}$. The conditional expectation in this model is calculated as $\mu_t = \frac{1}{p_t} - 1 = \exp(x_t' B)$. The probability is used in the calculation of the log-likelihood via the dgeom() function.

Several PMFs of the Geometric distribution are shown in the following figure:

alm() with distribution="dgeom" returns μ_t in mu and in fitted.values. The scale does not contain anything useful, because the Geometric distribution has a time varying variance, which depends on the probability. The prediction intervals are produces using qgeom() function. Similarly to the other count distribution, resid() method returns y_t − μ_t.

An example of application

Round up the response variable for the next example:

xreg[,1] <- round(abs(xreg[,1]))
inSample <- xreg[1:180,]
outSample <- xreg[-c(1:180),]

Negative Binomial distribution:

ourModel <- alm(y~x1+x2, data=inSample, distribution="dnbinom")
summary(ourModel)
#> Response variable: y
#> Distribution used in the estimation: Negative Binomial with size=33.8042
#> Loss function used in estimation: likelihood
#> Coefficients:
#>             Estimate Std. Error Lower 2.5% Upper 97.5%  
#> (Intercept)   5.9876     0.1417     5.7080      6.2671 *
#> x1           -0.0009     0.0046    -0.0100      0.0081  
#> x2            0.0034     0.0026    -0.0018      0.0086  
#> 
#> Error standard deviation: 82.4676
#> Sample size: 180
#> Number of estimated parameters: 4
#> Number of degrees of freedom: 176
#> Information criteria:
#>      AIC     AICc      BIC     BICc 
#> 2116.695 2116.924 2129.467 2130.060

And an example with predefined size:

ourModel <- alm(y~x1+x2, data=inSample, distribution="dnbinom", size=30)
summary(ourModel)
#> Response variable: y
#> Distribution used in the estimation: Negative Binomial with size=30
#> Loss function used in estimation: likelihood
#> Coefficients:
#>             Estimate Std. Error Lower 2.5% Upper 97.5%  
#> (Intercept)   6.0010     0.1501     5.7048      6.2972 *
#> x1           -0.0010     0.0048    -0.0105      0.0086  
#> x2            0.0032     0.0028    -0.0023      0.0087  
#> 
#> Error standard deviation: 82.4643
#> Sample size: 180
#> Number of estimated parameters: 4
#> Number of degrees of freedom: 176
#> Information criteria:
#>      AIC     AICc      BIC     BICc 
#> 2117.149 2117.378 2129.921 2130.514

Cumulative Logistic distribution

We have previously discussed the density function of logistic distribution. The standardised cumulative distribution function used in alm() is: where q̂_t = x_t′A is the conditional mean of the level, underlying the probability. This value is then used in the likelihood in order to estimate the parameters of the model. The error term of the model is calculated using the formula: This way the error varies from −∞ to ∞ and is equal to zero, when u_t = 0.5.

The alm() function with distribution="plogis" returns q_t in mu, standard deviation, calculated using the respective errors in scale and the probability p̂_t based on in fitted.values. resid() method returns the errors discussed above. predict() method produces point forecasts and the intervals for the probability of occurrence. The intervals use the assumption of normality of the error term, generating respective quantiles (based on the estimated q_t and variance of the error) and then transforming them into the scale of probability using Logistic CDF. This method for intervals calculation is approximate and should not be considered as a final solution!

Cumulative Normal distribution

The case of cumulative Normal distribution is quite similar to the cumulative Logistic one. The transformation is done using the standard normal CDF: where q_t = x_t′A. Similarly to the Logistic CDF, the estimated probability is used in the likelihood in order to estimate the parameters of the model. The error term is calculated using the standardised quantile function of Normal distribution: It acts similar to the error from the Logistic distribution, but is based on the different set of functions. Its CDF has similar shapes to the logit:

Similar to the Logistic CDF, the alm() function with distribution="pnorm" returns q_t in mu, standard deviation, calculated based on the errors in scale and the probability p̂_t based on in fitted.values. resid() method returns the errors discussed above. predict() method produces point forecasts and the intervals for the probability of occurrence. The intervals are also approximate and use the same principle as in Logistic CDF.