Package 'greybox'

Description

Function produces error measures for the provided object and the holdout values of the response variable. Note that instead of parameters x, test, the function accepts the vector of values in holdout. Also, the parameters d and D are not supported - MASE is always calculated via division by first differences.

Usage

## S3 method for class 'greybox'
accuracy(object, holdout = NULL, ...)

## S3 method for class 'predict.greybox'
accuracy(object, holdout = NULL, ...)

Arguments

Details

The function is a wrapper for the measures function and is implemented for convenience.

Author(s)

Ivan Svetunkov, [email protected]

Examples

xreg <- cbind(rlaplace(100,10,3),rnorm(100,50,5))
xreg <- cbind(100+0.5*xreg[,1]-0.75*xreg[,2]+rlaplace(100,0,3),xreg,rnorm(100,300,10))
colnames(xreg) <- c("y","x1","x2","Noise")

ourModel <- alm(y~x1+x2+trend, xreg, subset=c(1:80), distribution="dlaplace")
predict(ourModel,xreg[-c(1:80),]) |>
   accuracy(xreg[-c(1:80),"y"])

Description

This is a simple method that returns the values of the response variable of the model

Usage

actuals(object, all = TRUE, ...)

## Default S3 method:
actuals(object, all = TRUE, ...)

## S3 method for class 'lm'
actuals(object, all = TRUE, ...)

## S3 method for class 'alm'
actuals(object, all = TRUE, ...)

## S3 method for class 'predict.greybox'
actuals(object, all = TRUE, ...)

Arguments

Value

The vector of the response variable.

Author(s)

Ivan Svetunkov, [email protected]

Examples

xreg <- cbind(rnorm(100,10,3),rnorm(100,50,5))
xreg <- cbind(100+0.5*xreg[,1]-0.75*xreg[,2]+rnorm(100,0,3),xreg,rnorm(100,300,10))
colnames(xreg) <- c("y","x1","x2","Noise")

ourModel <- stepwise(xreg)

actuals(ourModel)

Description

This function extracts AICc / BICc from models. It can be applied to wide variety of models that use logLik() and nobs() methods (including the popular lm, forecast, smooth classes).

Usage

AICc(object, ...)

BICc(object, ...)

Arguments

Details

AICc was proposed by Nariaki Sugiura in 1978 and is used on small samples for the models with normally distributed residuals. BICc was derived in McQuarrie (1999) and is used in similar circumstances.

IMPORTANT NOTE: both of the criteria can only be used for univariate models (regression models, ARIMA, ETS etc) with normally distributed residuals! In case of multivariate models, both criteria need to be modified. See Bedrick & Tsai (1994) for details.

Value

This function returns numeric value.

Author(s)

Ivan Svetunkov, [email protected]

References

• Burnham Kenneth P. and Anderson David R. (2002). Model Selection and Multimodel Inference. A Practical Information-Theoretic Approach. Springer-Verlag New York. DOI: [10.1007/b97636](http://dx.doi.org/10.1007/b97636).

• McQuarrie, A. D. (1999). A small-sample correction for the Schwarz SIC model selection criterion. Statistics & Probability Letters, 44(1), 79–86. [10.1016/S0167-7152(98)00294-6](https://doi.org/10.1016/S0167-7152(98)00294-6).

• McQuarrie A.D., A small-sample correction for the Schwarz SIC model selection criterion, Statistics & Probability Letters 44 (1999) pp.79-86. doi:10.1016/S0167-7152(98)00294-6

• Sugiura Nariaki (1978) Further analysts of the data by Akaike's information criterion and the finite corrections, Communications in Statistics - Theory and Methods, 7:1, 13-26, doi:10.1080/03610927808827599

• Bedrick, E. J., & Tsai, C.-L. (1994). Model Selection for Multivariate Regression in Small Samples. Biometrics, 50(1), 226. doi:10.2307/2533213

See Also

AIC, BIC

Examples

xreg <- cbind(rnorm(100,10,3),rnorm(100,50,5))
xreg <- cbind(100+0.5*xreg[,1]-0.75*xreg[,2]+rnorm(100,0,3),xreg,rnorm(100,300,10))
colnames(xreg) <- c("y","x1","x2","Noise")

ourModel <- stepwise(xreg)

AICc(ourModel)
BICc(ourModel)

Description

The function applies several models on the provided time series and identifies what type of demand it is based on an information criterion.

Usage

aid(y, ic = c("AICc", "AIC", "BICc", "BIC"), level = 0.99,
  loss = "likelihood", ...)

Arguments

Details

In the first step, function creates inter-demand intervals and fits a model with LOWESS of it assuming Geometric distribution. The outliers from this model are treated as potential stock outs.

In the second step, the function creates explanatory variables based on LOWESS of the original data, then applies Normal, Normal + Bernoulli models and selects the one that has the lowest IC. Based on that, it decides what type of demand the data corresponds to: regular or intermittent. Finally, if the data is count, the function will identify that.

Value

Class "aid" is returned, which contains:

• y - The original data;

• models - All fitted models;

• ICs - Values of information criteria;

• type - The type of the identified demand;

• stockouts - List with start and end ids of potential stockouts;

• new - Binary showing whether the data start with the abnormal number of zeroes. Must be a new product then;

• obsolete - Binary showing whether the data ends with the abnormal number of zeroes. Must be product that was discontinued (obsolete).

Author(s)

Ivan Svetunkov, [email protected]

Examples

# Data from Poisson distribution
y <- rpois(120, 0.7)
aid(y)

Description

Function estimates model based on the selected distribution

Usage

alm(formula, data, subset, na.action, distribution = c("dnorm", "dlaplace",
  "ds", "dgnorm", "dlogis", "dt", "dalaplace", "dlnorm", "dllaplace", "dls",
  "dlgnorm", "dbcnorm", "dinvgauss", "dgamma", "dexp", "dfnorm", "drectnorm",
  "dpois", "dnbinom", "dbinom", "dgeom", "dbeta", "dlogitnorm", "plogis",
  "pnorm"), loss = c("likelihood", "MSE", "MAE", "HAM", "LASSO", "RIDGE",
  "ROLE"), occurrence = c("none", "plogis", "pnorm"), scale = NULL,
  orders = c(0, 0, 0), parameters = NULL, fast = FALSE, ...)

Arguments

Details

This is a function, similar to lm, but using likelihood for the cases of several non-normal distributions. These include:

1. dnorm - Normal distribution,

2. dlaplace - Laplace distribution,

3. ds - S-distribution,

4. dgnorm - Generalised Normal distribution,

5. dlogis - Logistic Distribution,

6. dt - T-distribution,

7. dalaplace - Asymmetric Laplace distribution,

8. dlnorm - Log-Normal distribution,

9. dllaplace - Log-Laplace distribution,

10. dls - Log-S distribution,

11. dlgnorm - Log-Generalised Normal distribution,

12. dfnorm - Folded normal distribution,

13. drectnorm - Rectified normal distribution,

14. dbcnorm - Box-Cox normal distribution,

15. dinvgauss - Inverse Gaussian distribution,

16. dgamma - Gamma distribution,

17. dexp - Exponential distribution,

18. dlogitnorm - Logit-normal distribution,

19. dbeta - Beta distribution,

20. dpois - Poisson Distribution,

21. dnbinom - Negative Binomial Distribution,

22. dbinom - Binomial Distribution,

23. dgeom - Geometric Distribution,

24. plogis - Cumulative Logistic Distribution,

25. pnorm - Cumulative Normal distribution.

This function can be considered as an analogue of glm, but with the focus on time series. This is why, for example, the function has orders parameter for ARIMA and produces time series analysis plots with plot(alm(...)).

This function is slower than lm, because it relies on likelihood estimation of parameters, hessian calculation and matrix multiplication. So think twice when using distribution="dnorm" here.

The estimation is done via the maximisation of likelihood of a selected distribution, so the number of estimated parameters always includes the scale. Thus the number of degrees of freedom of the model in case of alm will typically be lower than in the case of lm.

See more details and examples in the vignette for "ALM": vignette("alm","greybox")

Value

Function returns model - the final model of the class "alm", which contains:

• coefficients - estimated parameters of the model,

• FI - Fisher Information of parameters of the model. Returned only when FI=TRUE,

• fitted - fitted values,

• residuals - residuals of the model,

• mu - the estimated location parameter of the distribution,

• scale - the estimated scale parameter of the distribution. If a formula was provided for scale, then an object of class "scale" will be returned.

• distribution - distribution used in the estimation,

• logLik - log-likelihood of the model. Only returned, when loss="likelihood" or loss="ROLE" and in several special cases of distribution and loss combinations (e.g. loss="MSE", distribution="dnorm"),

• loss - the type of the loss function used in the estimation,

• lossFunction - the loss function, if the custom is provided by the user,

• lossValue - the value of the loss function,

• res - the output of the optimisation (nloptr function),

• df.residual - number of degrees of freedom of the residuals of the model,

• df - number of degrees of freedom of the model,

• call - how the model was called,

• rank - rank of the model,

• data - data used for the model construction,

• terms - terms of the data. Needed for some additional methods to work,

• occurrence - the occurrence model used in the estimation,

• B - the value of the optimised parameters. Typically, this is a duplicate of coefficients,

• other - the list of all the other parameters either passed to the function or estimated in the process, but not included in the standard output (e.g. alpha for Asymmetric Laplace),

• timeElapsed - the time elapsed for the estimation of the model.

Author(s)

Ivan Svetunkov, [email protected]

See Also

stepwise, lmCombine, xregTransformer

Examples

### An example with mtcars data and factors
mtcars2 <- within(mtcars, {
   vs <- factor(vs, labels = c("V", "S"))
   am <- factor(am, labels = c("automatic", "manual"))
   cyl  <- factor(cyl)
   gear <- factor(gear)
   carb <- factor(carb)
})
# The standard model with Log-Normal distribution
ourModel <- alm(mpg~., mtcars2[1:30,], distribution="dlnorm")
summary(ourModel)
plot(ourModel)
# Produce table based on the output for LaTeX
xtable(summary(ourModel))

# Produce predictions with the one sided interval (upper bound)
predict(ourModel, mtcars2[-c(1:30),], interval="p", side="u")

# Model with heteroscedasticity (scale changes with the change of qsec)
ourModel <- alm(mpg~., mtcars2[1:30,], scale=~qsec)

### Artificial data for the other examples
xreg <- cbind(rlaplace(100,10,3),rnorm(100,50,5))
xreg <- cbind(100+0.5*xreg[,1]-0.75*xreg[,2]+rlaplace(100,0,3),xreg,rnorm(100,300,10))
colnames(xreg) <- c("y","x1","x2","Noise")

# An example with Laplace distribution
ourModel <- alm(y~x1+x2+trend, xreg, subset=c(1:80), distribution="dlaplace")
summary(ourModel)
plot(predict(ourModel,xreg[-c(1:80),]))

# And another one with Asymmetric Laplace distribution (quantile regression)
# with optimised alpha
ourModel <- alm(y~x1+x2, xreg, subset=c(1:80), distribution="dalaplace")

# An example with AR(1) order
ourModel <- alm(y~x1+x2, xreg, subset=c(1:80), distribution="dnorm", orders=c(1,0,0))
summary(ourModel)
plot(predict(ourModel,xreg[-c(1:80),]))

### Examples with the count data
xreg[,1] <- round(exp(xreg[,1]-70),0)

# Negative Binomial distribution
ourModel <- alm(y~x1+x2, xreg, subset=c(1:80), distribution="dnbinom")
summary(ourModel)
predict(ourModel,xreg[-c(1:80),],interval="p",side="u")

# Poisson distribution
ourModel <- alm(y~x1+x2, xreg, subset=c(1:80), distribution="dpois")
summary(ourModel)
predict(ourModel,xreg[-c(1:80),],interval="p",side="u")


### Examples with binary response variable
xreg[,1] <- round(xreg[,1] / (1 + xreg[,1]),0)

# Logistic distribution (logit regression)
ourModel <- alm(y~x1+x2, xreg, subset=c(1:80), distribution="plogis")
summary(ourModel)
plot(predict(ourModel,xreg[-c(1:80),],interval="c"))

# Normal distribution (probit regression)
ourModel <- alm(y~x1+x2, xreg, subset=c(1:80), distribution="pnorm")
summary(ourModel)
plot(predict(ourModel,xreg[-c(1:80),],interval="p"))

Description

Function returns the matrix of measures of association for different types of variables.

Usage

association(x, y = NULL, use = c("na.or.complete", "complete.obs",
  "everything", "all.obs"), method = c("auto", "pearson", "spearman",
  "kendall", "cramer"))

assoc(x, y = NULL, use = c("na.or.complete", "complete.obs", "everything",
  "all.obs"), method = c("auto", "pearson", "spearman", "kendall", "cramer"))

Arguments

Details

The function looks at the types of the variables and calculates different measures depending on the result:

• If both variables are numeric, then Pearson's correlation is calculated;

• If both variables are categorical, then Cramer's V is calculated;

• Finally, if one of the variables is categorical, and the other is numeric, then multiple correlation is returned.

After that the measures are wrapped up in a matrix.

Function also calculates the p-values associated with the respective measures (see the return).

See details in the vignette "Marketing analytics with greybox": vignette("maUsingGreybox","greybox")

assoc() is just a short name for the association{}.

Value

The following list of values is returned:

• value - Matrix of the coefficients of association;

• p.value - The p-values for the parameters;

• type - The matrix of the types of measures of association.

Author(s)

Ivan Svetunkov, [email protected]

See Also

table, tableplot, spread, cramer, mcor

Examples

association(mtcars)

Description

These are the basic methods for the alm and greybox models that extract coefficients, their covariance matrix, confidence intervals or generating the summary of the model. If the non-likelihood related loss was used in the process, then it is recommended to use bootstrap (which is slow, but more reliable).

Usage

## S3 method for class 'greybox'
coef(object, bootstrap = FALSE, ...)

## S3 method for class 'alm'
confint(object, parm, level = 0.95, bootstrap = FALSE, ...)

## S3 method for class 'scale'
confint(object, parm, level = 0.95, bootstrap = FALSE, ...)

## S3 method for class 'alm'
vcov(object, bootstrap = FALSE, ...)

## S3 method for class 'scale'
vcov(object, bootstrap = FALSE, ...)

## S3 method for class 'alm'
summary(object, level = 0.95, bootstrap = FALSE, ...)

Arguments

Details

The coef() method returns the vector of parameters of the model. If bootstrap=TRUE, then the coefficients are calculated as the mean values of the bootstrapped ones.

The vcov() method returns the covariance matrix of parameters. If bootstrap=TRUE, then the bootstrap is done using coefbootstrap function

The confint() constructs the confidence intervals for parameters. Once again, this can be done using bootstrap=TRUE.

Finally, the summary() returns the table with parameters, their standard errors, confidence intervals and general information about the model.

Value

Depending on the used method, different values are returned.

Author(s)

Ivan Svetunkov, [email protected]

See Also

alm, coefbootstrap

Examples

# An example with ALM
ourModel <- alm(mpg~., mtcars, distribution="dlnorm")
coef(ourModel)
vcov(ourModel)
confint(ourModel)
summary(ourModel)

Description

The function does the bootstrap for parameters of models and returns covariance matrix together with the original bootstrapped data.

Usage

coefbootstrap(object, nsim = 1000, size = floor(0.75 * nobs(object)),
  replace = FALSE, prob = NULL, parallel = FALSE, method = c("dsr",
  "cr"), ...)

## S3 method for class 'lm'
coefbootstrap(object, nsim = 1000, size = floor(0.75 *
  nobs(object)), replace = FALSE, prob = NULL, parallel = FALSE,
  method = c("dsr", "cr"), ...)

## S3 method for class 'alm'
coefbootstrap(object, nsim = 1000, size = floor(0.75 *
  nobs(object)), replace = FALSE, prob = NULL, parallel = FALSE,
  method = c("dsr", "cr"), ...)

Arguments

Details

The function applies the same model as in the provided object on a smaller sample in order to get the estimates of parameters and capture the uncertainty about them. This is a simple implementation of the case resampling, which assumes that the observations are independent.

Value

Class "bootstrap" is returned, which contains:

• vcov - the covariance matrix of parameters;

• coefficients - the matrix with the bootstrapped coefficients.

• nsim - number of runs done;

• size - the sample size used in the bootstrap;

• replace - whether the sampling was done with replacement;

• prob - a vector of probability weights used in the process;

• parallel - whether the calculations were done in parallel;

• model - the name of the model used (the name of the function);

• timeElapsed - the time that was spend on the calculations.

Author(s)

Ivan Svetunkov, [email protected]

See Also

alm

Examples

# An example with ALM
ourModel <- alm(mpg~., mtcars, distribution="dlnorm", loss="HAM")
# A fast example with 10 iterations. Use at least 100 to get better results
coefbootstrap(ourModel, nsim=10)

Description

Function calculates Cramer's V for two categorical variables based on the table function

Usage

cramer(x, y, use = c("na.or.complete", "complete.obs", "everything",
  "all.obs"), unbiased = TRUE)

Arguments

Details

The function calculates Cramer's V and also returns the associated statistics from Chi-Squared test with the null hypothesis of independence of the two variables.

See details in the vignette "Marketing analytics with greybox": vignette("maUsingGreybox","greybox")

Value

The following list of values is returned:

• value - The value of Cramer's V;

• statistic - The value of Chi squared statistic associated with the Cramer's V;

• p.value - The p-value of Chi squared test associated with the Cramer's V;

• df - The number of degrees of freedom from the test.

Author(s)

Ivan Svetunkov, [email protected]

References

Wicher Bergsma (2013), A bias-correction for Cramér's V and Tschuprow's T. Journal of the Korean Statistical Society, 42, pp. 323-328. doi:10.1016/j.jkss.2012.10.002.

See Also

table, tableplot, spread, mcor, association

Examples

cramer(mtcars$am, mtcars$gear)

Description

Density, cumulative distribution, quantile functions and random number generation for the Asymmetric Laplace distribution with the location parameter mu, scale and the asymmetry parameter alpha.

Usage

dalaplace(q, mu = 0, scale = 1, alpha = 0.5, log = FALSE)

palaplace(q, mu = 0, scale = 1, alpha = 0.5)

qalaplace(p, mu = 0, scale = 1, alpha = 0.5)

ralaplace(n = 1, mu = 0, scale = 1, alpha = 0.5)

Arguments

Details

When mu=0 and scale=1, the Laplace distribution becomes standardized. The distribution has the following density function:

f(x) = alpha (1-alpha) / scale exp(-(x-mu)/scale (alpha - I(x<=mu))),

where I(.) is the indicator function (equal to 1 if the condition is satisfied and zero otherwise).

When alpha=0.5, then the distribution becomes Symmetric Laplace, where scale = 1/2 MAE.

This distribution function aligns with the quantile estimates of parameters (Geraci & Bottai, 2007).

Finally, both palaplace and qalaplace are returned for the lower tail of the distribution.

Value

Depending on the function, various things are returned (usually either vector or scalar):

• dalaplace returns the density function value for the provided parameters.

• palaplace returns the value of the cumulative function for the provided parameters.

• qalaplace returns quantiles of the distribution. Depending on what was provided in p, mu and scale, this can be either a vector or a matrix, or an array.

• ralaplace returns a vector of random variables generated from the Laplace distribution. Depending on what was provided in mu and scale, this can be either a vector or a matrix or an array.

Author(s)

Ivan Svetunkov, [email protected]

References

• Geraci Marco, Bottai Matteo (2007). Quantile regression for longitudinal data using the asymmetric Laplace distribution. Biostatistics (2007), 8, 1, pp. 140-154 doi:10.1093/biostatistics/kxj039

• Yu, K., & Zhang, J. (2005). A three-parameter asymmetric laplace distribution and its extension. Communications in Statistics - Theory and Methods, 34, 1867-1879. doi:10.1080/03610920500199018

See Also

Distributions

Examples

x <- dalaplace(c(-100:100)/10, 0, 1, 0.2)
plot(x, type="l")

x <- palaplace(c(-100:100)/10, 0, 1, 0.2)
plot(x, type="l")

qalaplace(c(0.025,0.975), 0, c(1,2), c(0.2,0.3))

x <- ralaplace(1000, 0, 1, 0.2)
hist(x)

Description

Density, cumulative distribution, quantile functions and random number generation for the distribution that becomes normal after the Box-Cox transformation. Note that this is based on the original Box-Cox paper.

Usage

dbcnorm(q, mu = 0, sigma = 1, lambda = 0, log = FALSE)

pbcnorm(q, mu = 0, sigma = 1, lambda = 0)

qbcnorm(p, mu = 0, sigma = 1, lambda = 0)

rbcnorm(n = 1, mu = 0, sigma = 1, lambda = 0)

Arguments

Details

The distribution has the following density function:

f(y) = y^(lambda-1) 1/sqrt(2 pi) exp(-((y^lambda-1)/lambda -mu)^2 / (2 sigma^2))

Both pbcnorm and qbcnorm are returned for the lower tail of the distribution.

In case of lambda=0, the values of the log normal distribution are returned. In case of lambda=1, the values of the normal distribution are returned with mu=mu+1.

All the functions are defined for non-negative values only.

Value

Depending on the function, various things are returned (usually either vector or scalar):

• dbcnorm returns the density function value for the provided parameters.

• pbcnorm returns the value of the cumulative function for the provided parameters.

• qbcnorm returns quantiles of the distribution. Depending on what was provided in p, mu and sigma, this can be either a vector or a matrix, or an array.

• rbcnorm returns a vector of random variables generated from the bcnorm distribution. Depending on what was provided in mu and sigma, this can be either a vector or a matrix or an array.

Author(s)

Ivan Svetunkov, [email protected]

References

• Box, G. E., & Cox, D. R. (1964). An Analysis of Transformations. Journal of the Royal Statistical Society. Series B (Methodological), 26(2), 211–252. Retrieved from https://www.jstor.org/stable/2984418

See Also

Distributions

Examples

x <- dbcnorm(c(-1000:1000)/200, 0, 1, 1)
plot(c(-1000:1000)/200, x, type="l")

x <- pbcnorm(c(-1000:1000)/200, 0, 1, 1)
plot(c(-1000:1000)/200, x, type="l")

qbcnorm(c(0.025,0.975), 0, c(1,2), 1)

x <- rbcnorm(1000, 0, 1, 1)
hist(x)

Description

Functions to detect, when Daylight Saving Time and leap year start and finish

Usage

detectdst(object)

detectleap(object)

Arguments

Details

The detectdst function detects, when the change for the DST starts and ends. It assumes that the frequency of data is not lower than hourly. The detectleap function does similar for the leap year, but flagging the 29th of February as a starting and to the 28th of February next year as the ending dates.

In order for the methods to work, the object needs to be of either zoo / xts or POSIXt class and should contain valid dates.

Value

List containing:

• start - data frame with id (number of observation) and the respective dates, when the DST / leap year start;

• end - data frame with id and dates, when DST / leap year end.

Author(s)

Ivan Svetunkov, [email protected]

See Also

xregExpander, temporaldummy, outlierdummy

Examples

# Generate matrix with monthly dummies for a zoo object
x <- as.POSIXct("2004-01-01")+0:(365*24*8)*60*60
detectdst(x)
detectleap(x)

Description

Function produces coefficients of determination for the provided data

Usage

determination(xreg, bruteforce = TRUE, ...)

determ(object, ...)

Arguments

Details

The function calculates coefficients of determination (aka R^2) between all the provided variables. The higher the coefficient for a variable is, the higher the potential multicollinearity effect in the model with the variable will be. Coefficients of determination are connected directly to Variance Inflation Factor (VIF): VIF = 1 / (1 - determination). Arguably it is easier to interpret, because it is restricted with (0, 1) bounds. The multicollinearity can be considered as serious, when determination > 0.9 (which corresponds to VIF > 10).

The method determ can be applied to wide variety of classes, including lm, glm and alm.

See details in the vignette "Marketing analytics with greybox": vignette("maUsingGreybox","greybox")

Value

Function returns the vector of determination coefficients.

Author(s)

Ivan Svetunkov, [email protected]

See Also

cor, mcor, stepwise

Examples

### Simple example
xreg <- cbind(rnorm(100,10,3),rnorm(100,50,5))
xreg <- cbind(100+0.5*xreg[,1]-0.75*xreg[,2]+rnorm(100,0,3),xreg,rnorm(100,300,10))
colnames(xreg) <- c("x1","x2","x3","Noise")
determination(xreg)

Description

Density, cumulative distribution, quantile functions and random number generation for the folded normal distribution with the location parameter mu and the scale sigma (which corresponds to standard deviation in normal distribution).

Usage

dfnorm(q, mu = 0, sigma = 1, log = FALSE)

pfnorm(q, mu = 0, sigma = 1)

qfnorm(p, mu = 0, sigma = 1)

rfnorm(n = 1, mu = 0, sigma = 1)

Arguments

Details

The distribution has the following density function:

f(x) = 1/sqrt(2 pi) (exp(-(x-mu)^2 / (2 sigma^2)) + exp(-(x+mu)^2 / (2 sigma^2)))

Both pfnorm and qfnorm are returned for the lower tail of the distribution.

Value

Depending on the function, various things are returned (usually either vector or scalar):

• dfnorm returns the density function value for the provided parameters.

• pfnorm returns the value of the cumulative function for the provided parameters.

• qfnorm returns quantiles of the distribution. Depending on what was provided in p, mu and sigma, this can be either a vector or a matrix, or an array.

• rfnorm returns a vector of random variables generated from the fnorm distribution. Depending on what was provided in mu and sigma, this can be either a vector or a matrix or an array.

Author(s)

Ivan Svetunkov, [email protected]

References

• Wikipedia page on folded normal distribution: https://en.wikipedia.org/wiki/Folded_normal_distribution.

See Also

Distributions

Examples

x <- dfnorm(c(-1000:1000)/200, 0, 1)
plot(x, type="l")

x <- pfnorm(c(-1000:1000)/200, 0, 1)
plot(x, type="l")

qfnorm(c(0.025,0.975), 0, c(1,2))

x <- rfnorm(1000, 0, 1)
hist(x)

Description

Density, cumulative distribution, quantile functions and random number generation for the Generalised Normal distribution with the location mu, a scale and a shape parameters.

Usage

dgnorm(q, mu = 0, scale = 1, shape = 1, log = FALSE)

pgnorm(q, mu = 0, scale = 1, shape = 1, lower.tail = TRUE,
  log.p = FALSE)

qgnorm(p, mu = 0, scale = 1, shape = 1, lower.tail = TRUE,
  log.p = FALSE)

rgnorm(n, mu = 0, scale = 1, shape = 1)

Arguments

Details

A generalized normal random variable $x$ with parameters location $\mu$, scale $s > 0$ and shape $\beta > 0$ has density:

$p(x) = \beta exp{-(|x - \mu|/s)^\beta}/(2s \Gamma(1/\beta)).$

The mean and variance of $x$ are $\mu$ and $s^2 \Gamma(3/\beta)/\Gamma(1/\beta)$, respectively.

The function are based on the functions from gnorm package of Maryclare Griffin (package has been abandoned since 2018).

The quantile and cumulative functions use uniform approximation for cases shape>100. This is needed, because otherwise it is not possible to calculate the values correctly due to scale^(shape)=Inf in R.

Author(s)

Maryclare Griffin and Ivan Svetunkov

Source

dgnorm, pgnorm, qgnorm andrgnorm are all parametrized as in Version 1 of the Generalized Normal Distribution Wikipedia page, which uses the parametrization given by in Nadarajah (2005). The same distribution was described much earlier by Subbotin (1923) and named the exponential power distribution by Box and Tiao (1973).

References

• Box, G. E. P. and G. C. Tiao. "Bayesian inference in Statistical Analysis." Addison-Wesley Pub. Co., Reading, Mass (1973).

• Nadarajah, Saralees. "A generalized normal distribution." Journal of Applied Statistics 32.7 (2005): 685-694.

• Subbotin, M. T. "On the Law of Frequency of Error." Matematicheskii Sbornik 31.2 (1923): 206-301.

See Also

Distributions

Examples

# Density function values for standard normal distribution
x <- dgnorm(seq(-1, 1, length.out = 100), 0, sqrt(2), 2)
plot(x, type="l")

#CDF of standard Laplace
x <- pgnorm(c(-100:100), 0, 1, 1)
plot(x, type="l")

# Quantiles of S distribution
qgnorm(c(0.025,0.975), 0, 1, 0.5)

# Random numbers from a distribution with shape=10000 (approximately uniform)
x <- rgnorm(1000, 0, 1, 1000)
hist(x)

Description

The greybox package implements several distribution functions. In this document, I list the main functions and provide links to related resources.

Details

• Generalised normal distribution (with a kurtosis parameter by Nadarajah, 2005): qgnorm, dgnorm, pgnorm, rgnorm.

• S distribution (a special case of Generalised Normal with shape=0.5): qs, ds, ps, rs.

• Laplace distribution (special case of Generalised Normal with shape=1): qlaplace, dlaplace, plaplace, rlaplace.

• Asymmetric Laplace distribution (Yu & Zhang, 2005): qalaplace, dalaplace, palaplace, ralaplace.

• Logit Normal distribution (Mead, 1965): qlogitnorm, dlogitnorm, plogitnorm, rlogitnorm.

• Box-Cox Normal distribution (Box & Cox, 1964): qbcnorm, dbcnorm, pbcnorm, rbcnorm.

• Folded Normal distribution: qfnorm, dfnorm, pfnorm, rfnorm.

• Rectified Normal distribution: qrectnorm, drectnorm, prectnorm, rrectnorm.

• Three parameter Log Normal distribution (Sangal & Biswas, 1970): qtplnorm, dtplnorm, ptplnorm, rtplnorm.

Author(s)

Ivan Svetunkov, [email protected]

Ivan Svetunkov, [email protected]

References

• Nadarajah, Saralees. "A generalized normal distribution." Journal of Applied Statistics 32.7 (2005): 685-694.

• Wikipedia page on Laplace distribution: https://en.wikipedia.org/wiki/Laplace_distribution.

• Yu, K., & Zhang, J. (2005). A three-parameter asymmetric laplace distribution and its extension. Communications in Statistics - Theory and Methods, 34, 1867-1879. doi:10.1080/03610920500199018

• Mead, R. (1965). A Generalised Logit-Normal Distribution. Biometrics, 21 (3), 721–732. doi: 10.2307/2528553

• Box, G. E., & Cox, D. R. (1964). An Analysis of Transformations. Journal of the Royal Statistical Society. Series B (Methodological), 26(2), 211–252. Retrieved from https://www.jstor.org/stable/2984418

• Sangal, B. P., & Biswas, A. K. (1970). The 3-Parameter Distribution Applications in Hydrology. Water Resources Research, 6(2), 505–515. doi:10.1029/WR006i002p00505

See Also

Distributions from the stats package.

Description

Density, cumulative distribution, quantile functions and random number generation for the Laplace distribution with the location parameter mu and the scale parameter (which is equal to Mean Absolute Error, aka Mean Absolute Deviation).

Usage

dlaplace(q, mu = 0, scale = 1, log = FALSE)

plaplace(q, mu = 0, scale = 1)

qlaplace(p, mu = 0, scale = 1)

rlaplace(n = 1, mu = 0, scale = 1)

Arguments

Details

When mu=0 and scale=1, the Laplace distribution becomes standardized. The distribution has the following density function:

f(x) = 1/(2 scale) exp(-abs(x-mu) / scale)

Both plaplace and qlaplace are returned for the lower tail of the distribution.

Value

Depending on the function, various things are returned (usually either vector or scalar):

• dlaplace returns the density function value for the provided parameters.

• plaplace returns the value of the cumulative function for the provided parameters.

• qlaplace returns quantiles of the distribution. Depending on what was provided in p, mu and scale, this can be either a vector or a matrix, or an array.

• rlaplace returns a vector of random variables generated from the Laplace distribution. Depending on what was provided in mu and scale, this can be either a vector or a matrix or an array.

Author(s)

Ivan Svetunkov, [email protected]

References

• Wikipedia page on Laplace distribution: https://en.wikipedia.org/wiki/Laplace_distribution.

See Also

Distributions

Examples

x <- dlaplace(c(-100:100)/10, 0, 1)
plot(x, type="l")

x <- plaplace(c(-100:100)/10, 0, 1)
plot(x, type="l")

qlaplace(c(0.025,0.975), 0, c(1,2))

x <- rlaplace(1000, 0, 1)
hist(x)

Description

Density, cumulative distribution, quantile functions and random number generation for the distribution that becomes normal after the Logit transformation.

Usage

dlogitnorm(q, mu = 0, sigma = 1, log = FALSE)

plogitnorm(q, mu = 0, sigma = 1)

qlogitnorm(p, mu = 0, sigma = 1)

rlogitnorm(n = 1, mu = 0, sigma = 1)

Arguments

Details

The distribution has the following density function:

f(y) = 1/(sqrt(2 pi) y (1-y)) exp(-(logit(y) -mu)^2 / (2 sigma^2))

where y is in (0, 1) and logit(y) =log(y/(1-y)).

Both plogitnorm and qlogitnorm are returned for the lower tail of the distribution.

All the functions are defined for the values between 0 and 1.

Value

Depending on the function, various things are returned (usually either vector or scalar):

• dlogitnorm returns the density function value for the provided parameters.

• plogitnorm returns the value of the cumulative function for the provided parameters.

• qlogitnorm returns quantiles of the distribution. Depending on what was provided in p, mu and sigma, this can be either a vector or a matrix, or an array.

• rlogitnorm returns a vector of random variables generated from the logitnorm distribution. Depending on what was provided in mu and sigma, this can be either a vector or a matrix or an array.

Author(s)

Ivan Svetunkov, [email protected]

References

• Mead, R. (1965). A Generalised Logit-Normal Distribution. Biometrics, 21 (3), 721–732. doi: 10.2307/2528553

See Also

Distributions

Examples

x <- dlogitnorm(c(-1000:1000)/200, 0, 1)
plot(c(-1000:1000)/200, x, type="l")

x <- plogitnorm(c(-1000:1000)/200, 0, 1)
plot(c(-1000:1000)/200, x, type="l")

qlogitnorm(c(0.025,0.975), 0, c(1,2))

x <- rlogitnorm(1000, 0, 1)
hist(x)

Description

Density, cumulative distribution, quantile functions and random number generation for the Rectified Normal distribution.

Usage

drectnorm(q, mu = 0, sigma = 1, log = FALSE)

prectnorm(q, mu = 0, sigma = 1)

qrectnorm(p, mu = 0, sigma = 1)

rrectnorm(n = 1, mu = 0, sigma = 1)

Arguments

Details

If x~N(mu, sigma^2) then y = max(0, x) follows Rectified Normal distribution: y~RectN(mu, sigma^2), which can be written as:

f_y = 1(x<=0) F_x(mu, sigma) + 1(x>0) f_x(x, mu, sigma),

where F_x is the cumulative distribution function and f_x is the probability density function of normal distribution.

Both prectnorm and qrectnorm are returned for the lower tail of the distribution.

All the functions are defined for non-negative values only.

Value

Depending on the function, various things are returned (usually either vector or scalar):

• drectnorm returns the density function value for the provided parameters.

• prectnorm returns the value of the cumulative function for the provided parameters.

• qrectnorm returns quantiles of the distribution. Depending on what was provided in p, mu and sigma, this can be either a vector or a matrix, or an array.

• rrectnorm returns a vector of random variables generated from the RectN distribution. Depending on what was provided in mu and sigma, this can be either a vector or a matrix or an array.

Author(s)

Ivan Svetunkov, [email protected]

See Also

Distributions

Examples

x <- drectnorm(c(-1000:1000)/200, 0, 1)
plot(c(-1000:1000)/200, x, type="l")

x <- prectnorm(c(-1000:1000)/200, 0, 1)
plot(c(-1000:1000)/200, x, type="l")

qrectnorm(c(0.025,0.975), 0, c(1,2))

x <- rrectnorm(1000, 0, 1)
hist(x)

Description

Density, cumulative distribution, quantile functions and random number generation for the S distribution with the location parameter mu and a scaling parameter scale.

Usage

ds(q, mu = 0, scale = 1, log = FALSE)

ps(q, mu = 0, scale = 1)

qs(p, mu = 0, scale = 1)

rs(n = 1, mu = 0, scale = 1)

Arguments

Details

When mu=0 and ham=2, the S distribution becomes standardized with scale=1 (this is because scale=ham/2). The distribution has the following density function:

f(x) = 1/(4 scale^2) exp(-sqrt(abs(x-mu)) / scale)

The S distribution has fat tails and large excess.

Both ps and qs are returned for the lower tail of the distribution.

Value

Depending on the function, various things are returned (usually either vector or scalar):

• ds returns the density function value for the provided parameters.

• ps returns the value of the cumulative function for the provided parameters.

• qs returns quantiles of the distribution. Depending on what was provided in p, mu and scale, this can be either a vector or a matrix, or an array.

• rs returns a vector of random variables generated from the S distribution. Depending on what was provided in mu and scale, this can be either a vector or a matrix or an array.

Author(s)

Ivan Svetunkov, [email protected]

See Also

Distributions

Examples

x <- ds(c(-1000:1000)/10, 0, 1)
plot(x, type="l")

x <- ps(c(-1000:1000)/10, 0, 1)
plot(x, type="l")

qs(c(0.025,0.975), 0, 1)

x <- rs(1000, 0, 1)
hist(x)

Description

The function implements a bootstrap inspired by the Maximum Entropy Bootstrap

Usage

dsrboot(y, nsim = 100, intermittent = FALSE, type = c("multiplicative",
  "additive"), kind = c("nonparametric", "parametric"), lag = frequency(y),
  sd = NULL, scale = TRUE)

## S3 method for class 'dsrboot'
plot(x, sorted = FALSE, legend = TRUE, ...)

Arguments

Details

The "Data Shape Replication" bootstrap reproduces the shape of the original time series by creating randomness around it. It is done in the following steps:

1. Sort the data in the ascending order, recording the original order of elements; 2. Take first differences of the original data and sort them; 3. Generate random numbers from the uniform distribution between 0 and 1; 4. Get the smoothed differences that correspond to the random numbers (randomly extract empirical quantiles). This way we take the empirical density into account when selecting the differences; 5. Add the random differences to the sorted series from (1) to get a new time series; 6. Sort the new time series in the ascending order; 7. Reorder (6) based on the initial order of series; 8. Centre the data around the original series; 9. Scale the data to make sure that the variance is constant over time.

If the multiplicative bootstrap is used then logarithms of the sorted series are used and at the very end, the exponent of the resulting data is taken. This way the discrepancies in the data have similar scale no matter what the level of the original series is. In case of the additive bootstrap, the trended series will be noisier when the level of series is low.

Value

The function returns:

• call - the call that was used originally;

• data - the original data used in the function;

• boot - the matrix with the new series in columns and observations in rows.

• type - type of the bootstrap used.

• sd - the value of sd used in case of parameteric bootstrap.

• scale - whether the scaling was needed.

• smooth - the smoothed ordered actual data.

Author(s)

Ivan Svetunkov, [email protected]

References

Vinod HD, López-de-Lacalle J (2009). "Maximum Entropy Bootstrap for Time Series: The meboot R Package." Journal of Statistical Software, 29(5), 1–19. doi:10.18637/jss.v029.i05.

Examples

dsrboot(AirPassengers) |> plot()

Description

Density, cumulative distribution, quantile functions and random number generation for the 3 parameter log normal distribution with the location parameter mu, scale sigma (which corresponds to standard deviation in normal distribution) and shifting parameter shift.

Usage

dtplnorm(q, mu = 0, sigma = 1, shift = 0, log = FALSE)

ptplnorm(q, mu = 0, sigma = 1, shift = 0)

qtplnorm(p, mu = 0, sigma = 1, shift = 0)

rtplnorm(n = 1, mu = 0, sigma = 1, shift = 0)

Arguments

Details

The distribution has the following density function:

f(x) = 1/(x-a) 1/sqrt(2 pi) exp(-(log(x-a)-mu)^2 / (2 sigma^2))

Both ptplnorm and qtplnorm are returned for the lower tail of the distribution.

The function is based on the lnorm functions from stats package, introducing the shift parameter.

Value

Depending on the function, various things are returned (usually either vector or scalar):

• dtplnorm returns the density function value for the provided parameters.

• ptplnorm returns the value of the cumulative function for the provided parameters.

• qtplnorm returns quantiles of the distribution. Depending on what was provided in p, mu and sigma, this can be either a vector or a matrix, or an array.

• rtplnorm returns a vector of random variables generated from the tplnorm distribution. Depending on what was provided in mu and sigma, this can be either a vector or a matrix or an array.

Author(s)

Ivan Svetunkov, [email protected]

References

• Sangal, B. P., & Biswas, A. K. (1970). The 3-Parameter Distribution Applications in Hydrology. Water Resources Research, 6(2), 505–515. doi:10.1029/WR006i002p00505

See Also

Distributions

Examples

x <- dtplnorm(c(-1000:1000)/200, 0, 1, 1)
plot(c(-1000:1000)/200, x, type="l")

x <- ptplnorm(c(-1000:1000)/200, 0, 1, 1)
plot(c(-1000:1000)/200, x, type="l")

qtplnorm(c(0.025,0.975), 0, c(1,2), 1)

x <- rtplnorm(1000, 0, 1, 1)
hist(x)

Description

This function allows extracting error type from any model.

Usage

errorType(object, ...)

Arguments

Details

errorType extracts the type of error from the model (either additive or multiplicative).

Value

Either "A" for additive error or "M" for multiplicative. All the other functions return strings of character.

Author(s)

Ivan Svetunkov, [email protected]

Examples

xreg <- cbind(rnorm(100,10,3),rnorm(100,50,5))
xreg <- cbind(100+0.5*xreg[,1]-0.75*xreg[,2]+rnorm(100,0,3),xreg,rnorm(100,300,10))
colnames(xreg) <- c("y","x1","x2","Noise")
ourModel <- alm(y~x1+x2,as.data.frame(xreg))

errorType(ourModel)

Description

Functions extract scale and the standard error of the residuals. Mainly needed for the work with the model estimated via sm.

Usage

extractScale(object, ...)

## Default S3 method:
extractScale(object, ...)

## S3 method for class 'greybox'
extractScale(object, ...)

extractSigma(object, ...)

## Default S3 method:
extractSigma(object, ...)

## S3 method for class 'greybox'
extractSigma(object, ...)

Arguments

Details

In case of a simpler model, the functions will return the scalar using sigma() method. If the scale model was estimated, extractScale() and extractSigma() will return the conditional scale and the conditional standard error of the residuals respectively.

Value

One of the two is returned, depending on the type of estimated model:

• Scalar from sigma() method if the variance is assumed to be constant.

• Vector of values if the scale model was estimated

Author(s)

Ivan Svetunkov, [email protected]

See Also

sm

Examples

# Generate the data
xreg <- cbind(rnorm(100,10,3),rnorm(100,50,5))
xreg <- cbind(100+0.5*xreg[,1]-0.75*xreg[,2]+sqrt(exp(0.8+0.2*xreg[,1]))*rnorm(100,0,1),
              xreg,rnorm(100,300,10))
colnames(xreg) <- c("y","x1","x2","Noise")

# Estimate the location and scale model
ourModel <- alm(y~., xreg, scale=~x1+x2)

# Extract scale
extractScale(ourModel)
# Extract standard error
extractSigma(ourModel)

Description

The function makes a standard linear graph using the provided actuals and forecasts.

Usage

graphmaker(actuals, forecast, fitted = NULL, lower = NULL, upper = NULL,
  level = NULL, legend = TRUE, cumulative = FALSE, vline = TRUE,
  parReset = TRUE, ...)

Arguments

Details

Function uses the provided data to construct a linear graph. It is strongly advised to use ts objects to define the start of each of the vectors. Otherwise the data may be plotted incorrectly.

Value

Function does not return anything.

Author(s)

Ivan Svetunkov

See Also

ts

Examples

xreg <- cbind(y=rnorm(100,100,10),x=rnorm(100,10,10))
almModel <- alm(y~x, xreg, subset=c(1:90))
values <- predict(almModel, newdata=xreg[-c(1:90),], interval="prediction")

graphmaker(xreg[,1],values$mean,fitted(values))
graphmaker(xreg[,1],values$mean,fitted(values),legend=FALSE)
graphmaker(xreg[,1],values$mean,fitted(values),legend=FALSE,lower=values$lower,upper=values$upper)

# Produce the necessary ts objects from an arbitrary vectors
actuals <- ts(c(1:10), start=c(2000,1), frequency=4)
forecast <- ts(c(11:15),start=end(actuals)[1]+end(actuals)[2]*deltat(actuals),
               frequency=frequency(actuals))
graphmaker(actuals,forecast)

# This should work as well
graphmaker(c(1:10),c(11:15))

# This way you can add additional elements to the plot
graphmaker(c(1:10),c(11:15), parReset=FALSE)
points(c(1:15))
# But don't forget to do dev.off() in order to reset the plotting area afterwards

Description

Toolbox for working with univariate models for purposes of analysis and forecasting

Details

The following functions are included in the package:

• AICc and BICc - AIC / BIC corrected for the sample size.

• pointLik - point likelihood of the function.

• pAIC, pAICc, pBIC, pBICc - point versions of respective information criteria.

• coefbootstrap - Method that uses a simple implementation of the case resampling to get bootstrapped estimates of parameters of the model.

• dsrboot - Bootstrap inspired by the meboot package, that creates bootstraped series based on the provided one.

• determination - Coefficients of determination between different exogenous variables.

• temporaldummy - Matrix with seasonal dummy variables.

• outlierdummy - Matrix with dummies for outliers.

• alm - Advanced Linear Model - regression, estimated using likelihood with specified distribution (e.g. Laplace or Chi-Squared).

• sm - Scale Model - Regression model for scale of distribution (e.g. for Variance of Normal distribution). Requires an lm() or alm() model.

• stepwise - Stepwise based on information criteria and partial correlations. Efficient and fast.

• xregExpander - Function that expands the provided data into the data with lags and leads.

• xregTransformer - Function produces mathematical transformations of the variables, such as taking logarithms, square roots etc.

• xregMultiplier - Function produces cross-products of the matrix of the provided variables.

• lmCombine - Function combines lm models from the estimated based on information criteria weights.

• lmDynamic - Dynamic regression based on point AIC.

• ro - Rolling origin evaluation.

• Distributions - document explaining the distribution functions of the greybox package.

• spread - function that produces scatterplots / boxplots / tableplots, depending on the types of variables.

• assoc - function that calculates measures of association, depending on the types of variables.

Author(s)

Ivan Svetunkov, [email protected]

See Also

stepwise, lmCombine

Examples

xreg <- cbind(rnorm(100,10,3),rnorm(100,50,5))
xreg <- cbind(100+0.5*xreg[,1]-0.75*xreg[,2]+rnorm(100,0,3),xreg,rnorm(100,300,10))
colnames(xreg) <- c("y","x1","x2","Noise")

stepwise(xreg)

Description

hm() function estimates half moment from some predefined constant C. ham() estimates the Half Absolute Moment. asymmetry() function returns Asymmetry coefficient, while extremity() returns the coefficient of Extremity, both based on hm(). Finally, cextremity() returns the Complex Extremity coefficient, based on hm().

Usage

hm(x, C = mean(x, na.rm = TRUE), ...)

ham(x, C = mean(x, na.rm = TRUE), ...)

asymmetry(x, C = mean(x, na.rm = TRUE), ...)

extremity(x, C = mean(x, na.rm = TRUE), ...)

cextremity(x, C = mean(x, na.rm = TRUE), ...)

Arguments

Details

NA values of x are excluded on the first step of calculation.

Value

A complex variable is returned for the hm() and cextremity() functions, and real values are returned for ham(), asymmetry() and extremity().

Author(s)

Ivan Svetunkov, [email protected]

References

• Svetunkov I., Kourentzes N., Svetunkov S. "Half Central Moment for Data Analysis". Working Paper of Department of Management Science, Lancaster University, 2023:3, 1–21.

Examples

x <- rnorm(100,0,1)
hm(x)
ham(x)
asymmetry(x)
extremity(x)
cextremity(x)

Description

The function implants the scale model into the location model. It currently works with alm / adam and sm() method.

Usage

implant(location, scale, ...)

Arguments

Details

The function is needed in order to treat the scale of model correctly in the methods like forecast().

Value

The model of the same class as the location model, but with scale from the estimated model via sm(). This is needed to produce appropriate forecasts in case of scale model and to take into account the correct number of estimated parameters.

Author(s)

Ivan Svetunkov, [email protected]

See Also

alm, adam, sm

Examples

xreg <- cbind(rnorm(100,10,3),rnorm(100,50,5))
xreg <- cbind(100+0.5*xreg[,1]-0.75*xreg[,2]+sqrt(exp(0.8+0.2*xreg[,1]))*rnorm(100,0,1),
              xreg,rnorm(100,300,10))
colnames(xreg) <- c("y","x1","x2","Noise")

# Estimate the location model
ourModel <- alm(y~.,xreg)
# Estimate the scale model
ourScale <- sm(ourModel,formula=~x1+x2)
# Implant scale into location
ourModel <- implant(ourModel, ourScale)

Description

Functions to check if an object is of the specified class

Usage

is.greybox(x)

is.alm(x)

is.occurrence(x)

is.greyboxC(x)

is.greyboxD(x)

is.rollingOrigin(x)

is.rmc(x)

is.scale(x)

Arguments

Details

The list of functions includes:

• is.greybox() tests if the object was produced by a greybox function (e.g. alm / stepwise / lmCombine / lmDynamic);

• is.alm() tests if the object was produced by alm() function;

• is.occurrence() tests if an occurrence part of the model was produced;

• is.greyboxC() tests if the object was produced by lmCombine() function;

• is.greyboxD() tests if the object was produced by lmDynamic() function;

• is.rmc() tests if the object was produced by rmc() function;

• is.rollingOrigin() tests if the object was produced by ro() function;

• is.scale() tests if the object is of the class "scale" (produced by alm or sm in case of heteroscedastic model);

Value

TRUE if this is the specified class and FALSE otherwise.

Author(s)

Ivan Svetunkov, [email protected]

Examples

xreg <- cbind(rlaplace(100,10,3),rnorm(100,50,5))
xreg <- cbind(100+0.5*xreg[,1]-0.75*xreg[,2]+rlaplace(100,0,3),xreg,rnorm(100,300,10))
colnames(xreg) <- c("y","x1","x2","Noise")

ourModel <- alm(y~x1+x2, xreg, distribution="dnorm")

is.alm(ourModel)
is.greybox(ourModel)
is.greyboxC(ourModel)
is.greyboxD(ourModel)

Description

Function combines parameters of linear regressions of the first variable on all the other provided data.

Usage

lmCombine(data, ic = c("AICc", "AIC", "BIC", "BICc"), bruteforce = FALSE,
  silent = TRUE, formula = NULL, subset = NULL,
  distribution = c("dnorm", "dlaplace", "ds", "dgnorm", "dlogis", "dt",
  "dalaplace", "dlnorm", "dllaplace", "dls", "dlgnorm", "dbcnorm", "dinvgauss",
  "dgamma", "dexp", "dfnorm", "drectnorm", "dpois", "dnbinom", "dbeta",
  "dlogitnorm", "plogis", "pnorm"), parallel = FALSE, ...)

Arguments

Details

The algorithm uses alm() to fit different models and then combines the models based on the selected IC. The parameters are combined so that if they are not present in some of models, it is assumed that they are equal to zero. Thus, there is a shrinkage effect in the combination.

Some details and examples of application are also given in the vignette "Greybox": vignette("greybox","greybox")

Value

Function returns model - the final model of the class "greyboxC". The list of variables:

• coefficients - combined parameters of the model,

• vcov - combined covariance matrix of the model,

• fitted - the fitted values,

• residuals - residual of the model,

• distribution - distribution used in the estimation,

• logLik - combined log-likelihood of the model,

• IC - the values of the combined information criterion,

• ICType - the type of information criterion used,

• df.residual - number of degrees of freedom of the residuals of the combined model,

• df - number of degrees of freedom of the combined model,

• importance - importance of the parameters,

• combination - the table, indicating which variables were used in every model construction and what were the weights for each model,

• timeElapsed - the time elapsed for the estimation of the model.

Author(s)

Ivan Svetunkov, [email protected]

References

• Burnham Kenneth P. and Anderson David R. (2002). Model Selection and Multimodel Inference. A Practical Information-Theoretic Approach. Springer-Verlag New York. DOI: [10.1007/b97636](http://dx.doi.org/10.1007/b97636).

• McQuarrie, A. D. (1999). A small-sample correction for the Schwarz SIC model selection criterion. Statistics & Probability Letters, 44(1), 79–86. [10.1016/S0167-7152(98)00294-6](https://doi.org/10.1016/S0167-7152(98)00294-6).

See Also

step, xregExpander, stepwise

Examples

### Simple example
xreg <- cbind(rnorm(100,10,3),rnorm(100,50,5))
xreg <- cbind(100+0.5*xreg[,1]-0.75*xreg[,2]+rnorm(100,0,3),xreg,rnorm(100,300,10))
colnames(xreg) <- c("y","x1","x2","Noise")
inSample <- xreg[1:80,]
outSample <- xreg[-c(1:80),]
# Combine all the possible models
ourModel <- lmCombine(inSample,bruteforce=TRUE)
predict(ourModel,outSample)
plot(predict(ourModel,outSample))

### Fat regression example
xreg <- matrix(rnorm(5000,10,3),50,100)
xreg <- cbind(100+0.5*xreg[,1]-0.75*xreg[,2]+rnorm(50,0,3),xreg,rnorm(50,300,10))
colnames(xreg) <- c("y",paste0("x",c(1:100)),"Noise")
inSample <- xreg[1:40,]
outSample <- xreg[-c(1:40),]
# Combine only the models close to the optimal
ourModel <- lmCombine(inSample, ic="BICc",bruteforce=FALSE)
summary(ourModel)
plot(predict(ourModel, outSample))

# Combine in parallel - should increase speed in case of big data
## Not run: ourModel <- lmCombine(inSample, ic="BICc", bruteforce=TRUE, parallel=TRUE)
summary(ourModel)
plot(predict(ourModel, outSample))
## End(Not run)

Description

Function combines parameters of linear regressions of the first variable on all the other provided data using pAIC weights. This is an extension of the lmCombine function, which relies upon the idea that the combination weights might change over time.

Usage

lmDynamic(data, ic = c("AICc", "AIC", "BIC", "BICc"), bruteforce = FALSE,
  silent = TRUE, formula = NULL, subset = NULL,
  distribution = c("dnorm", "dlaplace", "ds", "dgnorm", "dlogis", "dt",
  "dalaplace", "dlnorm", "dllaplace", "dls", "dlgnorm", "dbcnorm", "dfnorm",
  "dinvgauss", "dgamma", "dpois", "dnbinom", "dlogitnorm", "plogis", "pnorm"),
  parallel = FALSE, lowess = TRUE, f = NULL, ...)

Arguments

Details

The algorithm uses alm() to fit different models and then combines the models based on the selected point IC. The combination weights are calculated for each observation based on the point IC and then smoothed via LOWESS if the respective parameter (lowess) is set to TRUE.

Some details and examples of application are also given in the vignette "Greybox": vignette("greybox","greybox")

Value

Function returns model - the final model of the class "greyboxD", which includes time varying parameters and dynamic importance of each variable. The list of variables:

• coefficients - the mean (over time) parameters of the model,

• vcov - the combined covariance matrix of the model,

• fitted - the fitted values,

• residuals - the residuals of the model,

• distribution - the distribution used in the estimation,

• logLik - the mean (over time) log-likelihood of the model,

• IC - dynamic values of the information criterion (pIC),

• ICType - the type of information criterion used,

• df.residual - mean number of degrees of freedom of the residuals of the model,

• df - mean number of degrees of freedom of the model,

• importance - dynamic importance of the parameters,

• call - call used in the function,

• rank - rank of the combined model,

• data - the data used in the model,

• mu - the location value of the distribution,

• scale - the scale parameter if alm() was used,

• coefficientsDynamic - table with parameters of the model, varying over the time,

• df.residualDynamic - dynamic df.residual,

• dfDynamic - dynamic df,

• weights - the dynamic weights for each model under consideration,

• timeElapsed - the time elapsed for the estimation of the model.

Author(s)

Ivan Svetunkov, [email protected]

References

• Burnham Kenneth P. and Anderson David R. (2002). Model Selection and Multimodel Inference. A Practical Information-Theoretic Approach. Springer-Verlag New York. DOI: [10.1007/b97636](http://dx.doi.org/10.1007/b97636).

• McQuarrie, A. D. (1999). A small-sample correction for the Schwarz SIC model selection criterion. Statistics & Probability Letters, 44(1), 79–86. [10.1016/S0167-7152(98)00294-6](https://doi.org/10.1016/S0167-7152(98)00294-6).

See Also

stepwise, lmCombine

Examples

### Simple example
xreg <- cbind(rnorm(100,10,3),rnorm(100,50,5))
xreg <- cbind(100+0.5*xreg[,1]-0.75*xreg[,2]+rnorm(100,0,3),xreg,rnorm(100,300,10))
colnames(xreg) <- c("y","x1","x2","Noise")
inSample <- xreg[1:80,]
outSample <- xreg[-c(1:80),]
# Combine all the possible models
ourModel <- lmDynamic(inSample,bruteforce=TRUE)
predict(ourModel,outSample)
plot(predict(ourModel,outSample))

Description

Function calculates multiple correlation between y and x, constructing a linear regression model

Usage

mcor(x, y, use = c("na.or.complete", "complete.obs", "everything",
  "all.obs"))

Arguments

Details

This is based on the linear regression model with the set of variables in x. The returned value is just a coefficient of multiple correlation from regression, the F-statistics of the model (thus testing the null hypothesis that all the parameters are equal to zero), the associated p-value and the degrees of freedom.

See details in the vignette "Marketing analytics with greybox": vignette("maUsingGreybox","greybox")

Value

The following list of values is returned:

• value - The value of the coefficient;

• statistic - The value of F-statistics associated with the parameter;

• p.value - The p-value of F-statistics associated with the parameter;

• df.residual - The number of degrees of freedom for the residuals;

• df - The number of degrees of freedom for the data.

Author(s)

Ivan Svetunkov, [email protected]

See Also

table, tableplot, spread, cramer, association

Examples

mcor(mtcars$am, mtcars$mpg)

Description

Functions allow to calculate different types of errors for point and interval predictions:

1. ME - Mean Error,

2. MAE - Mean Absolute Error,

3. MSE - Mean Squared Error,

4. MRE - Mean Root Error (Kourentzes, 2014),

5. MIS - Mean Interval Score (Gneiting & Raftery, 2007),

6. MPE - Mean Percentage Error,

7. MAPE - Mean Absolute Percentage Error (See Svetunkov, 2017 for the critique),

8. MASE - Mean Absolute Scaled Error (Hyndman & Koehler, 2006),

9. RMSSE - Root Mean Squared Scaled Error (used in M5 Competition),

10. rMAE - Relative Mean Absolute Error (Davydenko & Fildes, 2013),

11. rRMSE - Relative Root Mean Squared Error,

12. rAME - Relative Absolute Mean Error,

13. rMIS - Relative Mean Interval Score,

14. sMSE - Scaled Mean Squared Error (Petropoulos & Kourentzes, 2015),

15. sPIS- Scaled Periods-In-Stock (Wallstrom & Segerstedt, 2010),

16. sCE - Scaled Cumulative Error,

17. sMIS - Scaled Mean Interval Score,

18. GMRAE - Geometric Mean Relative Absolute Error.

Usage

ME(holdout, forecast, na.rm = TRUE)

MAE(holdout, forecast, na.rm = TRUE)

MSE(holdout, forecast, na.rm = TRUE)

MRE(holdout, forecast, na.rm = TRUE)

MIS(holdout, lower, upper, level = 0.95, na.rm = TRUE)

MPE(holdout, forecast, na.rm = TRUE)

MAPE(holdout, forecast, na.rm = TRUE)

MASE(holdout, forecast, scale, na.rm = TRUE)

RMSSE(holdout, forecast, scale, na.rm = TRUE)

rMAE(holdout, forecast, benchmark, na.rm = TRUE)

rRMSE(holdout, forecast, benchmark, na.rm = TRUE)

rAME(holdout, forecast, benchmark, na.rm = TRUE)

rMIS(holdout, lower, upper, benchmarkLower, benchmarkUpper, level = 0.95,
  na.rm = TRUE)

sMSE(holdout, forecast, scale, na.rm = TRUE)

sPIS(holdout, forecast, scale, na.rm = TRUE)

sCE(holdout, forecast, scale, na.rm = TRUE)

sMIS(holdout, lower, upper, scale, level = 0.95, na.rm = TRUE)

GMRAE(holdout, forecast, benchmark, na.rm = TRUE)

Arguments

Details

In case of sMSE, scale needs to be a squared value. Typical one – squared mean value of in-sample actuals.

If all the measures are needed, then measures function can help.

There are several other measures, see details of pinball and hm.

Value

All the functions return the scalar value.

Author(s)

Ivan Svetunkov, [email protected]

References

• Kourentzes N. (2014). The Bias Coefficient: a new metric for forecast bias https://kourentzes.com/forecasting/2014/12/17/the-bias-coefficient-a-new-metric-for-forecast-bias/

• Svetunkov, I. (2017). Naughty APEs and the quest for the holy grail. https://openforecast.org/2017/07/29/naughty-apes-and-the-quest-for-the-holy-grail/

• Fildes R. (1992). The evaluation of extrapolative forecasting methods. International Journal of Forecasting, 8, pp.81-98.

• Hyndman R.J., Koehler A.B. (2006). Another look at measures of forecast accuracy. International Journal of Forecasting, 22, pp.679-688.

• Petropoulos F., Kourentzes N. (2015). Forecast combinations for intermittent demand. Journal of the Operational Research Society, 66, pp.914-924.

• Wallstrom P., Segerstedt A. (2010). Evaluation of forecasting error measurements and techniques for intermittent demand. International Journal of Production Economics, 128, pp.625-636.

• Davydenko, A., Fildes, R. (2013). Measuring Forecasting Accuracy: The Case Of Judgmental Adjustments To Sku-Level Demand Forecasts. International Journal of Forecasting, 29(3), 510-522. doi:10.1016/j.ijforecast.2012.09.002

• Gneiting, T., & Raftery, A. E. (2007). Strictly proper scoring rules, prediction, and estimation. Journal of the American Statistical Association, 102(477), 359–378. doi:10.1198/016214506000001437

See Also

pinball, hm, measures

Examples

y <- rnorm(100,10,2)
testForecast <- rep(mean(y[1:90]),10)

MAE(y[91:100],testForecast)
MSE(y[91:100],testForecast)

MPE(y[91:100],testForecast)
MAPE(y[91:100],testForecast)

# Measures from Petropoulos & Kourentzes (2015)
MASE(y[91:100],testForecast,mean(abs(y[1:90])))
sMSE(y[91:100],testForecast,mean(abs(y[1:90]))^2)
sPIS(y[91:100],testForecast,mean(abs(y[1:90])))
sCE(y[91:100],testForecast,mean(abs(y[1:90])))

# Original MASE from Hyndman & Koehler (2006)
MASE(y[91:100],testForecast,mean(abs(diff(y[1:90]))))

testForecast2 <- rep(y[91],10)
# Relative measures, from and inspired by Davydenko & Fildes (2013)
rMAE(y[91:100],testForecast2,testForecast)
rRMSE(y[91:100],testForecast2,testForecast)
rAME(y[91:100],testForecast2,testForecast)
GMRAE(y[91:100],testForecast2,testForecast)

#### Measures for the prediction intervals
# An example with mtcars data
ourModel <- alm(mpg~., mtcars[1:30,], distribution="dnorm")
ourBenchmark <- alm(mpg~1, mtcars[1:30,], distribution="dnorm")

# Produce predictions with the interval
ourForecast <- predict(ourModel, mtcars[-c(1:30),], interval="p")
ourBenchmarkForecast <- predict(ourBenchmark, mtcars[-c(1:30),], interval="p")

MIS(mtcars$mpg[-c(1:30)],ourForecast$lower,ourForecast$upper,0.95)
sMIS(mtcars$mpg[-c(1:30)],ourForecast$lower,ourForecast$upper,mean(mtcars$mpg[1:30]),0.95)
rMIS(mtcars$mpg[-c(1:30)],ourForecast$lower,ourForecast$upper,
       ourBenchmarkForecast$lower,ourBenchmarkForecast$upper,0.95)

### Also, see pinball function for other measures for the intervals

Description

Function calculates several error measures using the provided forecasts and the data for the holdout sample.

Usage

measures(holdout, forecast, actual, digits = NULL, benchmark = c("naive",
  "mean"))

Arguments

Value

The functions returns the named vector of errors:

• ME,

• MAE,

• MSE

• MPE,

• MAPE,

• MASE,

• sMAE,

• RMSSE,

• sMSE,

• sCE,

• rMAE,

• rRMSE,

• rAME,

• asymmetry,

• sPIS.

For the details on these errors, see Errors.

Author(s)

Ivan Svetunkov, [email protected]

References

• Svetunkov, I. (2017). Naughty APEs and the quest for the holy grail. https://openforecast.org/2017/07/29/naughty-apes-and-the-quest-for-the-holy-grail/

• Fildes R. (1992). The evaluation of extrapolative forecasting methods. International Journal of Forecasting, 8, pp.81-98.

• Hyndman R.J., Koehler A.B. (2006). Another look at measures of forecast accuracy. International Journal of Forecasting, 22, pp.679-688.

• Petropoulos F., Kourentzes N. (2015). Forecast combinations for intermittent demand. Journal of the Operational Research Society, 66, pp.914-924.

• Wallstrom P., Segerstedt A. (2010). Evaluation of forecasting error measurements and techniques for intermittent demand. International Journal of Production Economics, 128, pp.625-636.

• Davydenko, A., Fildes, R. (2013). Measuring Forecasting Accuracy: The Case Of Judgmental Adjustments To Sku-Level Demand Forecasts. International Journal of Forecasting, 29(3), 510-522. doi:10.1016/j.ijforecast.2012.09.002

Examples

y <- rnorm(100,10,2)
ourForecast <- rep(mean(y[1:90]),10)

measures(y[91:100],ourForecast,y[1:90],digits=5)

Description

nparam() returns the number of estimated parameters in the model, while nvariate() returns number of variates for the response variable.

Usage

nparam(object, ...)

nvariate(object, ...)

Arguments

Details

nparam() is a very basic and a simple function which does what it says: extracts number of estimated parameters in the model. nvariate() returns number of variates (dimensions, columns) for the response variable (1 for the univariate regression).

Value

Both functions return numeric values.

Author(s)

Ivan Svetunkov, [email protected]

See Also

nobs, logLik

Examples

### Simple example
xreg <- cbind(rnorm(100,10,3),rnorm(100,50,5))
xreg <- cbind(100+0.5*xreg[,1]-0.75*xreg[,2]+rnorm(100,0,3),xreg,rnorm(100,300,10))
colnames(xreg) <- c("y","x1","x2","Noise")
ourModel <- lm(y~.,data=as.data.frame(xreg))

nparam(ourModel)
nvariate(ourModel)

Description

Function detects outliers and creates a matrix with dummy variables. Only point outliers are considered (no level shifts).

Usage

outlierdummy(object, ...)

## Default S3 method:
outlierdummy(object, level = 0.999, type = c("rstandard",
  "rstudent"), ...)

## S3 method for class 'alm'
outlierdummy(object, level = 0.999, type = c("rstandard",
  "rstudent"), ...)

Arguments

Details

The detection is done based on the type of distribution used and confidence level specified by user.

Value

The class "outlierdummy", which contains the list:

• outliers - the matrix with the dummy variables, flagging outliers;

• statistic - the value of the statistic for the normalised variable;

• id - the ids of the outliers (which observations have them);

• level - the confidence level used in the process;

• type - the type of the residuals used;

• errors - the errors used in the detection. In case of count distributions, probabilities are returned.

Author(s)

Ivan Svetunkov, [email protected]

See Also

influence.measures

Examples

# Generate the data with S distribution
xreg <- cbind(rnorm(100,10,3),rnorm(100,50,5))
xreg <- cbind(100+0.5*xreg[,1]-0.75*xreg[,2]+rs(100,0,3),xreg)
colnames(xreg) <- c("y","x1","x2")

# Fit the normal distribution model
ourModel <- alm(y~x1+x2, xreg, distribution="dnorm")

# Detect outliers
xregOutlierDummy <- outlierdummy(ourModel)

Description

This function returns a vector of AIC values for the in-sample observations

Usage

pAIC(object, ...)

pAICc(object, ...)

pBIC(object, ...)

pBICc(object, ...)

Arguments

Details

This is based on pointLik function. The formula for this is: pAIC_t = 2 * k - 2 * T * l_t , where k is the number of parameters, T is the number of observations and l_t is the point likelihood. This way we preserve the property that AIC = mean(pAIC).

Value

The function returns the vector of point AIC values.

Author(s)

Ivan Svetunkov, [email protected]

See Also

pointLik

Examples

xreg <- cbind(rnorm(100,10,3),rnorm(100,50,5))
xreg <- cbind(100+0.5*xreg[,1]-0.75*xreg[,2]+rnorm(100,0,3),xreg,rnorm(100,300,10))
colnames(xreg) <- c("y","x1","x2","Noise")
ourModel <- alm(y~x1+x2,as.data.frame(xreg))

pAICValues <- pAIC(ourModel)

mean(pAICValues)
AIC(ourModel)

Description

Function calculates partial correlations between the provided variables

Usage

pcor(x, y = NULL, use = c("na.or.complete", "complete.obs", "everything",
  "all.obs"), method = c("pearson", "spearman", "kendall"))

Arguments

Details

The calculation is done based on multiple linear regressions. The function calculates them for each pair of variables based on the residuals of linear models of those variables from the other variables in the dataset.

Value

The following list of values is returned:

• value - Matrix of the coefficients of partial correlations;

• p.value - The p-values for the parameters;

• method - The method used in the calculations.

Author(s)

Ivan Svetunkov, [email protected]