Simulation RStudio

Viewed 212

I'm using package fGarch in R to simulate a skewed normal distribution. The statistics I have are as follows:

Mean = 155

Sd = 35

Median = 150

I used the formula: (3*(mean-median))/sd to calculate skewness which came out to 0.427857 I used the following code:

x <- rsnorm(10000,mean = 155,sd=35,xi=0.427857)
plot(density(x))

However, the plot showed was negatively skewed even though my skewness parameter was positive. I ran the code a few times to ensure it wasn't a one-time thing. I can't seem to find the error (since I expected it to be positively skewed). Please help.

1 Answers

There are many definitions for skewness, such as the one you used Pearson's second skewness. The skewness Xi in fGarch is based on a paper (Fernandez & Steel 2000), and is an inverse scale factor in the positive and the negative directions.

As shown in the figure below, inverting skewness Xi produces the mirror image around zero. When Xi = 1, the distribution is symmetrical. When 0 < Xi < 1, negative skew. When Xi > 1, positive skew.

enter image description here

In your example, the Mean is skewed to the right of the Median. So we're looking for a positive skew (Xi > 1). By trial and error, the Xi could be located around 1.094095 . The median value given by the quantile function qsnorm(0.5, Mean, Sd, xi=1.094095) is indeed 150.

The simulated median is however a bit biased, which results in an unexpected Pearson's second skewness.

Mean <- 155
Sd <- 35
Median <- 150
3 * (Mean - Median) / Sd #Pearson's second skewness
#0.4285714

library(fGarch)
set.seed(1)
x <- rsnorm(100000, Mean, Sd, xi=1.094095)

mean(x)
#154.9353
sd(x)
#35.07564
median(x) #a bit higher than expected
#153.9029
3 * (mean(x) - median(x)) / sd(x) #biased from expected
#0.08830281

Finally, the code for locating Xi using General-purpose Optimization optim method.

fn <- function(log_Xi, Mean, Sd, Median) {
  abs(qsnorm(0.5, Mean, Sd, exp(log_Xi)) - Median)
  }

tem <- optim(0, fn, gr=NULL, Mean, Sd, Median)
Xi <- exp(tem$par)
#1.094095

And the code for the density plot above.

x_ <- seq(-5,5, by=0.01)
plot(x_, dsnorm(x_, xi=3), type='l', ylab='density', col='red', xlab='x', main='Skewness')
lines(x_, dsnorm(x_, xi=1), col='black', lty=2)
lines(x_, dsnorm(x_, xi=1/3), col='blue')
legend('topright', legend=c('3','1','1/3'), fill=c('red','black','blue'), title='Xi')
Related