3
$\begingroup$

I would like to simulate a multiple linear regression model using R.

If I have the skewness and kurtosis for the residuals, how can I do that?

$\endgroup$
  • $\begingroup$ Specifying skewness and kurtosis is not sufficient to determine the distribution. Can you explain more? $\endgroup$ – Glen_b May 21 '13 at 5:37
  • 2
    $\begingroup$ Furthermore, can you show what you've tried? $\endgroup$ – John May 21 '13 at 6:11
4
$\begingroup$

Although skewness and kurtosis do not determine the distribution of the residuals, you can choose a distribution (within mathematical constraints) that has zero mean, unit variance, and achieves the desired skewness and kurtosis (and then later rescale and shift it if desired). Perhaps the best known family of such distributions is the Pearson system. If such a distribution would be a suitable model for your residuals, then the simulation can proceed by:

  1. Choosing values of all explanatory variables $x_1, \ldots, x_p$ and forming them by columns into a matrix $X$.

  2. Choosing the coefficients $\beta_1, \ldots, \beta_p$.

  3. Generating iid variates $\varepsilon$ from the Pearson distribution with the given skewness and kurtosis (and mean and variance).

  4. Computing the dependent variable as $y = X\beta + \varepsilon$.

It is not guaranteed that the residuals from the model $y \sim X$ will have the given moments! The properties of the residuals depend on the method of fitting, among other things. But if the residuals are not terribly skew and are generated from a zero-mean distribution, then their distribution will generally cleave to the desired one.

Here is an R implementation and illustration. It uses the PearsonDS library to generate the errors.

library(PearsonDS)
set.seed(17)
n <- 320
beta <- c(-1,2,3)
moments <- c(0, 1, -1, 5)
x <- cbind(rep(1, n), matrix(rnorm((length(beta)-1)*n), nrow=n))
e <- rpearson(n, moments=moments)
y <- x %*% beta + e
fit <- lm(y ~ x - 1) # The constant is included in `x` already
par(mfrow=c(1,3))
hist(e, main="Errors")
hist(residuals(fit), main="Residuals")
plot(residuals(fit), e, main="Residuals vs. Errors")

Figures

The visual similarity of the histograms indicates that the randomly generated errors and the OLS residuals have qualitatively the same distribution (with a slightly negative skewness of $-1$ and slightly positive excess kurtosis of $5-3=2$). The tight scatterplot confirms this.

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.