ALM

ALM — Augmented Linear Model

Overview

The Augmented Linear Model (alm in R, ALM in Python) is the core estimator in greybox. It extends the standard linear model by supporting 26 distributions and 7 loss functions, enabling regression for a wide variety of data types — continuous, count, binary, bounded, and zero-inflated.

In R, alm() follows R's formula-based interface: alm(y ~ x1 + x2, data, distribution="dnorm").

In Python, ALM follows the scikit-learn estimator pattern: construct with ALM(distribution="dnorm"), then call model.fit(X, y). The formula module converts R-style formulas to design matrices.

Usage

R

library(greybox)

# Fit model
model <- alm(mpg ~ ., mtcars, distribution="dlnorm", subset=c(1:30))
summary(model)

# Predict
predict(model, newdata=tail(mtcars, 2), interval="prediction", level=c(0.9,0.95))

Python

import numpy as np
from greybox import ALM, formula, mtcars

# Parse formula and fit model
y_vec, X_mat = formula("mpg ~ .", mtcars.head(30))
model = ALM(distribution="dlnorm")
model.fit(X_mat, y_vec)
print(model.summary())

# Predict
_, X_new = formula("mpg ~ .", mtcars.tail(2))
model.predict(X_new, interval="prediction", level=[0.9,0.95])

Parameters

Parameters

Parameter	R	Python	Type	Default	Description
Formula	formula	—	formula	—	R formula, e.g. y ~ x1 + x2; in Python use formula() to get X, y
Data	data	—	data.frame	—	Data containing the variables; in Python pass X, y directly
Subset	subset	—	vector	NULL	Subset of observations to use
NA action	na.action	—	function	—	How to handle NAs
Distribution	distribution	distribution	str	"dnorm"	Distribution (see Distributions table)
Loss	loss	loss	str	"likelihood"	Loss function (see Loss Functions table)
Occurrence	occurrence	occurrence	str	"none"	Occurrence model: "none", "plogis", "pnorm"
Scale formula	scale	scale_formula	formula / array-like	NULL / None	Scale model formula
ARIMA orders	orders	orders	vector / tuple	c(0,0,0) / (0,0,0)	ARIMA orders (p, d, q); MA (q) not implemented, must be 0
Initial params	parameters	—	vector	NULL	Initial parameter values
Fast mode	fast	—	logical	FALSE	OLS fast estimation; not implemented in Python
Alpha	alpha (via ...)	alpha	float	NULL / None	Asymmetric Laplace quantile level (0–1)
Shape	shape (via ...)	shape	float	NULL / None	Generalized Normal shape parameter
Box-Cox lambda	lambdaBC (via ...)	lambda_bc	float	NULL / None	Box-Cox lambda (0–1)
Size	size (via ...)	size	float	NULL / None	Negative Binomial / Binomial size
Nu	nu (via ...)	nu	float	NULL / None	Student's t / Chi-squared degrees of freedom
Trim	—	trim	float	0.0	Trim proportion for ROLE loss
L1 lambda	—	lambda_l1	float	None	L1 regularization weight (LASSO)
L2 lambda	—	lambda_l2	float	None	L2 regularization weight (RIDGE)
Optimizer	—	nlopt_kargs	dict	{}	Optimizer settings (see Optimizer section)
Verbosity	—	verbose	int	0	Verbosity level

Distributions

All 26 distributions supported by ALM:

Continuous (location-scale)

Code	Name	Link	Extra Param	Use Case
dnorm	Normal	identity	—	General continuous data
dlaplace	Laplace	identity	—	Heavy-tailed, robust to outliers
ds	S (half-Laplace)	identity	—	Light-tailed data
dgnorm	Generalized Normal	identity	shape	Flexible tail weight
dlogis	Logistic	identity	—	Heavy tails, symmetric
dt	Student's t	identity	nu	Heavy tails, small samples
dalaplace	Asymmetric Laplace	identity	alpha	Quantile regression

Log-transformed continuous

Code	Name	Link	Extra Param	Use Case
dlnorm	Log-Normal	identity (log scale)	—	Positive, right-skewed data
dllaplace	Log-Laplace	identity (log scale)	—	Positive, heavy-tailed
dls	Log-S	identity (log scale)	—	Positive, light-tailed
dlgnorm	Log-Generalized Normal	identity (log scale)	shape	Positive, flexible tails

Bounded / special continuous

Code	Name	Link	Extra Param	Use Case
dfnorm	Folded Normal	identity	—	Absolute/non-negative data
drectnorm	Rectified Normal	identity	—	Zero-inflated non-negative
dbcnorm	Box-Cox Normal	identity	lambda_bc	Non-normal, power transform
dlogitnorm	Logit-Normal	identity (logit scale)	—	Data in (0, 1)
dbeta	Beta	identity	—	Data in (0, 1), two-part model

Non-negative continuous (log-link)

Code	Name	Link	Extra Param	Use Case
dinvgauss	Inverse Gaussian	log	—	Positive, right-skewed
dgamma	Gamma	log	—	Positive, right-skewed
dexp	Exponential	log	—	Time between events

Discrete (log-link)

Code	Name	Link	Extra Param	Use Case
dpois	Poisson	log	—	Count data (mean ≈ variance)
dnbinom	Negative Binomial	log	size	Overdispersed count data
dbinom	Binomial	log	—	Binary count data
dgeom	Geometric	log	—	Trials until first success
dchisq	Chi-squared	squared	nu	Sum of squared normals

Binary / occurrence

Code	Name	Link	Extra Param	Use Case
plogis	Logistic CDF	logistic	—	Binary classification
pnorm	Normal CDF (probit)	probit	—	Binary classification

Loss Functions

Loss	Description
likelihood	Maximum likelihood estimation (default)
MSE	Mean Squared Error
MAE	Mean Absolute Error
HAM	Half Absolute Moment
LASSO	Likelihood + L1 penalty (lambda_l1)
RIDGE	Likelihood + L2 penalty (lambda_l2)
ROLE	Robust Likelihood Estimation (trimmed, trim)

Attributes / Return Values

Attributes

Attribute	R	Python	Type	Description
Coefficients	coefficients	coefficients / coef	ndarray	Estimated coefficients
Intercept	—	intercept_	float	Intercept value
Fitted values	fitted.values / fitted	fitted	ndarray	Fitted values
Residuals	residuals	residuals	ndarray	Model residuals
Actuals	—	actuals / data	ndarray	Original response values
Location (mu)	mu	—	ndarray	Location parameter values (identity-link: same as fitted)
Scale	scale	scale	float	Scale parameter
Sigma	—	sigma	float	Residual standard error
Log-likelihood	logLik	loglik / log_lik	float	Log-likelihood
AIC	—	aic	float	Akaike Information Criterion
AICc	—	aicc	float	Corrected AIC
BIC	—	bic	float	Bayesian Information Criterion
BICc	—	bicc	float	Corrected BIC
Loss value	lossValue	loss_value	float	Final loss function value
Time elapsed	timeElapsed	time_elapsed	float	Computation time (seconds)
Residual df	df.residual	df_residual_	int	Residual degrees of freedom
Model df	df	—	int	Model degrees of freedom
Distribution	distribution	distribution_	str	Distribution name
Loss type	loss	loss_	str	Loss function name
Observations	—	nobs	int	Number of observations
Parameters	—	nparam / n_param	int / dict	Number of parameters; n_param returns {"number": nparam, "df": df_residual}
Formula	—	formula	str	Formula string (if provided)
Extra params	other	other_	dict	Extra parameter values (alpha, shape, etc.)
Fisher info	FI	—	matrix	Fisher Information matrix
ARIMA polynomial	—	arima_polynomial_	dict or None	AR+I polynomial coefficients keyed by lag name; None if no ARIMA
ARIMA string	—	arima_string_	str or None	Human-readable ARIMA order string, e.g. "ARIMA(1,0,0)"; None if no ARIMA

Methods

Method	R	Python	Description
summary	summary(model)	model.summary(level=0.95)	Summary with coefficients, std errors, p-values, significance
predict	predict(model, newdata, interval, level, side)	model.predict(X, interval, level, side)	Predictions with optional intervals
confint	confint(model, parm, level)	model.confint(parm, level)	Confidence intervals for parameters
vcov	vcov(model)	model.vcov()	Variance-covariance matrix
coef	coef(model)	model.coef (property)	Extract coefficients
sigma	sigma(model)	model.sigma (property)	Residual standard error
logLik	logLik(model)	model.loglik (property)	Log-likelihood
nobs	nobs(model)	model.nobs (property)	Number of observations
nparam	nparam(model)	model.nparam (property)	Number of parameters
score	—	model.score(X, y, metric)	Score using metric: "likelihood", "MSE", "MAE", "R2"
get_params	—	model.get_params()	Get constructor parameters (sklearn compat)
set_params	—	model.set_params(**params)	Set constructor parameters (sklearn compat)

PredictionResult (Python)

The predict() method returns a PredictionResult object.

Slots

Slot	Type	Description
mean	ndarray	Point forecasts
lower	ndarray or None	Lower interval bounds
upper	ndarray or None	Upper interval bounds
level	float or list	Confidence level(s) used
variances	ndarray or None	Variance estimates per observation
side	str	"both", "upper", or "lower"
interval	str	"none", "confidence", or "prediction"

Properties and methods

Method / Property	Description
to_dataframe()	Convert to pandas DataFrame
columns	Column names of DataFrame representation
shape	Shape of DataFrame representation
index	Index of DataFrame representation
values	Values of DataFrame representation
len(pred)	Number of observations
pred["mean"]	Column access by name

In R, predict() returns a list with $mean, $lower, $upper, $level, $variances.

Examples

Basic Normal Regression

# R
model <- alm(mpg ~ wt + hp, mtcars, distribution="dnorm")
summary(model)

# Python
import numpy as np
from greybox import ALM, formula

data = {"y": np.array([21,22,23,20,18,15,14,24,23,20]),
        "x1": np.array([2.6,2.8,2.3,3.2,3.4,3.8,4.0,2.2,2.5,3.0]),
        "x2": np.array([110,110,93,150,150,245,245,62,95,123])}

y, X = formula("y ~ x1 + x2", data)
model = ALM(distribution="dnorm")
model.fit(X, y)
print(model.summary())

Laplace Distribution (robust)

# R — robust to outliers
model <- alm(y ~ x, data, distribution="dlaplace")

# Python
model = ALM(distribution="dlaplace")
model.fit(X, y)

Prediction Intervals

# R — prediction interval
pred <- predict(model, newdata=new_data, interval="prediction", level=0.95)
pred$mean
pred$lower
pred$upper

# One-sided upper interval
pred <- predict(model, newdata=new_data, interval="prediction", level=0.95, side="upper")

# Python — prediction interval
pred = model.predict(X_new, interval="prediction", level=0.95)
print(pred.to_dataframe())

# One-sided upper interval
pred = model.predict(X_new, interval="prediction", level=0.95, side="upper")

# Multiple confidence levels
pred = model.predict(X_new, interval="prediction", level=[0.80, 0.95])

Note: Prediction and confidence interval bounds are computed using distribution-specific quantile functions to match R's predict.alm(). For example, dlnorm intervals are asymmetric and always positive; dpois/dnbinom intervals are integer-valued; plogis/pnorm bounds lie in (0, 1). Symmetric t-distribution bounds are used only for dnorm and dt.

Summary Interpretation

# Python
summary = model.summary()
print(summary)
# Output:
# Response variable: y
# Distribution used in the estimation: Normal
# Loss function type: likelihood
# Coefficients:
#                Estimate Std. Error Lower 2.5% Upper 97.5%
#  (Intercept)    5.0123     1.2345     2.5926      7.4320 *
#           x1    2.0045     0.1234     1.7626      2.2464 *
#
# Error standard deviation: 3.0123
# ...

ARIMA (Time Series) Models

ALM supports autoregressive and integrated (differencing) components via the orders=(p, d, q) parameter. Currently p (AR) and d (I) are supported; q (MA) must be 0. Lagged response values are appended automatically to the design matrix during fitting and prediction.

AR(1) model

# R
ts_data <- c(2.1, 3.4, 3.0, 4.2, 5.1, 4.8, 5.5, 6.0, 5.8, 6.5,
             7.1, 6.8, 7.5, 8.0, 7.8)
model_ar <- alm(y ~ 1, data.frame(y=ts_data), distribution="dnorm", orders=c(1,0,0))
summary(model_ar)
predict(model_ar, h=3, interval="prediction", level=0.9)

# Python
import numpy as np
from greybox import ALM
from greybox.formula import formula

ts = np.array([2.1, 3.4, 3.0, 4.2, 5.1, 4.8, 5.5, 6.0, 5.8, 6.5,
               7.1, 6.8, 7.5, 8.0, 7.8])
y, X = formula("y ~ 1", {"y": ts})  # intercept-only design matrix

model_ar = ALM(distribution="dnorm", orders=(1, 0, 0))
model_ar.fit(X, y)

print(model_ar.summary())
# Response variable: y
# Distribution used in the estimation: Normal
# Loss function type: likelihood
# ARIMA(1,0,0)
# Coefficients:
#                Estimate Std. Error Lower 2.5% Upper 97.5%
#  (Intercept)    1.xxxx     x.xxxx     x.xxxx      x.xxxx *
#        yLag1   -x.xxxx     x.xxxx    -x.xxxx     -x.xxxx *
# ...

print(model_ar.arima_string_)        # 'ARIMA(1,0,0)'
print(model_ar.arima_polynomial_)    # {'yLag1': -0.xx}

# Forecast 3 steps ahead with 90% prediction intervals
X_new = np.ones((3, 1))  # intercept column only (lags appended automatically)
pred = model_ar.predict(X_new, interval="prediction", level=0.9)
print(pred.to_dataframe())
#     mean    lower    upper
# 0   x.xx    x.xx     x.xx
# 1   x.xx    x.xx     x.xx
# 2   x.xx    x.xx     x.xx

ARI(1,1) model with exogenous regressor

# Python
import numpy as np
from greybox import ALM
from greybox.formula import formula

np.random.seed(42)
n = 40
x = np.linspace(0, 4, n)
y_raw = 1.5 * x + np.cumsum(np.random.randn(n) * 0.5) + 10.0
y2, X2 = formula("y ~ x", {"y": y_raw, "x": x})

model_ari = ALM(distribution="dnorm", orders=(1, 1, 0))
model_ari.fit(X2, y2)

print(model_ari.summary())
print(model_ari.arima_string_)      # 'ARIMA(1,1,0)'
print(model_ari.arima_polynomial_)  # {'yLag1': ..., 'yLag2': ...}

# Out-of-sample forecast (supply new x values; y placeholder not used)
x_new = np.array([4.1, 4.2, 4.3])
_, X_new2 = formula("y ~ x", {"y": np.zeros(3), "x": x_new})
pred2 = model_ari.predict(X_new2, interval="prediction", level=0.95)
print(pred2.to_dataframe())

ARDL Distributed Lag Models

An Autoregressive Distributed Lag (ARDL) model combines AR dynamics (orders parameter) with distributed lags of exogenous regressors created with the backshift operator B(). This allows modelling the full dynamic response of y to a one-unit shock in x over time.

B() — Backshift Operator

B(x, k) creates a lagged (k > 0) or lead (k < 0) copy of variable x. It is designed to be used directly inside formula strings.

Signature

# R
B(x, k, ...)          # k > 0: lag, k < 0: lead, k = 0: identity

# Python
from greybox import B
B(x, k, gaps="auto")  # same semantics

Parameters

Parameter	R	Python	Type	Default	Description
x	x	x	array-like	—	Input variable
k	k	k	int	—	Lag order; positive = lag, negative = lead, 0 = identity
gaps	... (via xregExpander)	gaps	str	"auto"	Boundary fill: "auto", "NAs", "zero", "naive", "extrapolate"

Simple lag/lead examples

# R
y <- rnorm(10)
B(y, 1)   # lag-1: B(y,1)[t] = y[t-1]
B(y, -1)  # lead-1: B(y,-1)[t] = y[t+1]
B(y, 2)   # lag-2

# Python
import numpy as np
from greybox import B

y = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
B(y, 1)   # array([1., 1., 2., 3., 4.])  (boundary extrapolated)
B(y, -1)  # array([2., 3., 4., 5., 5.])
B(y, 2)   # array([1., 1., 1., 2., 3.])

multipliers() — Dynamic Multipliers

multipliers() computes the impulse-response of y to a one-unit shock in x over forecast horizon h. It reads the distributed-lag coefficients (x, B(x,1), B(x,2), ...) from a fitted ALM model and combines them with the AR polynomial via the recurrence:

m_0 = β_0
m_s = β_s + Σ_{i=1}^{min(p,s)} φ_i · m_{s-i}

Signature

# R
multipliers(object, parm, h = 10)

# Python
from greybox import multipliers
multipliers(model, parm, h=10)

Parameters

Parameter	R	Python	Type	Default	Description
object / model	object	model	ALM	—	Fitted ALM model
parm	parm	parm	str	—	Variable name as it appears in the design matrix
h	h	h	int	10	Forecast horizon

Returns

Named numeric vector (R) or dict (Python): {"h1": m1, "h2": m2, ..., "hh": mh}.

Example: ARDL with distributed lags and AR(1)

# R — Seatbelts dataset (UK road casualties)
library(greybox)
model <- alm(drivers ~ kms + B(kms, 1) + B(kms, 2),
             Seatbelts, distribution = "dnorm", orders = c(1, 0, 0))
summary(model)

# Dynamic multipliers of kms over 10 periods
multipliers(model, "kms", h = 10)
# h1        h2        h3  ...
# 0.xxxx    0.xxxx    0.xxxx

# Python — synthetic ARDL(1,0,0) with distributed lags
import numpy as np
from greybox import ALM, formula, B, multipliers

np.random.seed(42)
n = 60
x = np.random.randn(n).cumsum()
y = 0.5 * x + 0.3 * B(x, 1) + 0.2 * B(x, 2) + np.random.randn(n) * 0.5
data = {"y": y, "x": x}

y_vec, X_mat = formula("y ~ x + B(x, 1) + B(x, 2)", data)
model = ALM(distribution="dnorm", orders=(1, 0, 0))
model.fit(X_mat, y_vec)
print(model.summary())
# Coefficients:
#                Estimate  Std. Error ...
#  (Intercept)  ...
#  x            ~0.50
#  B(x,1)       ~0.30
#  B(x,2)       ~0.20
#  yLag1        ...

# Dynamic multipliers: cumulative effect of a 1-unit shock to x
dm = multipliers(model, "x", h=10)
print(dm)
# {'h1': 0.50, 'h2': 0.77, 'h3': 0.82, ...}

Optimizer Configuration

ALM uses nlopt for optimization (nloptr in case of R). Settings can be customized via the nlopt_kargs parameter (Python) or ... arguments (R).

Python nlopt_kargs options

Key	Type	Default	Description
algorithm	str	"NLOPT_LN_NELDERMEAD"	Optimization algorithm
maxeval	int	40 * n_params	Maximum function evaluations
maxtime	float	600	Maximum time in seconds
xtol_rel	float	1e-6	Relative tolerance on parameters
xtol_abs	float	1e-8	Absolute tolerance on parameters
ftol_rel	float	1e-4	Relative tolerance on function value
ftol_abs	float	0	Absolute tolerance on function value
print_level	int	0	Verbosity (0=none, 3=full)

Available algorithms

Algorithm	Description
NLOPT_LN_NELDERMEAD	Nelder-Mead simplex (default, derivative-free)
NLOPT_LN_SBPLX	Subplex (improved Nelder-Mead)
NLOPT_LN_COBYLA	COBYLA (constrained)
NLOPT_LN_BOBYQA	BOBYQA (quadratic approx)
NLOPT_LN_NEWUOA	NEWUOA (unconstrained quadratic)
NLOPT_LN_PRAXIS	PRAXIS (gradient-free)
NLOPT_LD_LBFGS	L-BFGS (gradient-based)
NLOPT_LD_SLSQP	SLSQP (gradient-based, constrained)
NLOPT_LD_MMA	MMA (gradient-based, constrained)

# Example: custom optimizer settings
model = ALM(
    distribution="dnorm",
    nlopt_kargs={
        "algorithm": "NLOPT_LN_SBPLX",
        "maxeval": 1000,
        "xtol_rel": 1e-8,
        "print_level": 1,
    }
)

R vs Python Differences

Aspect	R	Python
Interface	alm(formula, data)	ALM().fit(X, y)
Formula	Built-in	Via formula() → X, y
Coefficients	coef(model) or model$coefficients	model.coef property
Intercept	Included in coef()	Separate model.intercept_
Predict	predict(model, newdata=df)	model.predict(X)
Summary	summary(model)	model.summary()
Optimizer	nloptr	nlopt
fast param	Available (TRUE for OLS)	Not implemented
sklearn compat	—	get_params(), set_params(), score()

References

• Svetunkov, I. (2023). Statistics for Business Analytics. https://openforecast.org/sba/
• Svetunkov, I. (2022). Forecasting and Analytics with the Augmented Dynamic Adaptive Model (ADAM). https://openforecast.org/adam/