**Edit after @eric_kernfeld answer.**
I'd like to do
1. Generate a time-series, for example, from a uniform distribution.
2. Transform of non-normal variable to standard normal distribution.
3. Fit an arima model to standard normal variable.
4. Simulate from the arima model with the fitted parameters (in this case errors should be standard normal).
5. Apply the back transformation that converts simulated arima output to marginal variable.
On the steps 2 and 5 I'm going to use the Johnson transformation and the back Johnson transformation respectively.
I have a time series and I'd like to do a simulation of log-returns using the normalization with [Johnson distribution][1].
library(JohnsonDistribution)
library(moments)
rm(list=ls(all=TRUE))
set.seed(1)
n <- 252
log_mydata1 <- runif(n, min=-0.06, max= 0.05)
m1 <- mean(log_mydata1)
var1 <- sd(log_mydata1)
sk1 <- skewness(log_mydata1)
k1 <- kurtosis(log_mydata1)
# Fitting Johnson distribution parameters and type
FitJohnsonDistribution(m1, var1, sk1, k1)
iType <- FitJohnsonDistribution(m1, var1, sk1, k1)[1]
gamma <- FitJohnsonDistribution(m1, var1, sk1, k1)[2]
delta <- FitJohnsonDistribution(m1, var1, sk1, k1)[3]
lambda <- FitJohnsonDistribution(m1, var1, sk1, k1)[4]
xi <- FitJohnsonDistribution(m1, var1, sk1, k1)[5]
# Applying Johnson transformation
z1 <- zJohnsonDistribution(log_mydata1, iType, gamma, delta, lambda , xi)
shapiro.test(z1) # W = 0.99374, p-value = 0.377
I have used [the code][2] and fitted my random data and log-returns were fitted with `ARIMA(3,0,2)` model. Then I applied some test to check the model quality.
# Fitting ARMA model
ArimaModelFit <- function(z)
{
final.aic <- Inf
final.order <- c(0,0,0)
for (p in 0:3)
for (q in 0:3)
{
if ( p == 0 && q == 0) {
next
}
arimaFit = tryCatch( arima(z, order=c(p, 0, q)),
error=function( err ) FALSE,
warning=function( err ) FALSE )
if( !is.logical( arimaFit ) ) {
current.aic <- AIC(arimaFit)
if (current.aic < final.aic) {
final.aic <- current.aic
final.order <- c(p, 0, q)
final.arima <- arima(z, order=final.order)
}
} else {
next
}
}
result <- list(aic=final.aic, order=final.order, arima=final.arima)
return(result)
} # function
f1 <- ArimaModelFit(z1)
rf1 <- residuals(f1$arima); shapiro.test(rf1) # W = 0.9944, p-value = 0.4785
Then I simulated data with the ARIMA model
#ARMA Simulation
sim <- arima.sim(list(order = c(3,0,2),
ar = c(f1$arima$coef[1], f1$arima$coef[2], f1$arima$coef[3]),
ma = c(f1$arima$coef[4], f1$arima$coef[5])), n = n)
Finally, I'd like to back the fitted data to marginal variable. I have applied the Kolmogorov-Smirnov test and plotted Cumulative Distribution Functions (CDFs) of marginal (`log-returns`) and simulated (`y`) data to check the quality of 2, 3, 4, 5 step of simulation (normalization, arima, back transformation).
# Applying inverse of Johnson transformation
y <- yJohnsonDistribution(sim, iType, gamma, delta, lambda , xi)
# Two-sample Kolmogorov-Smirnov test
ks.test(log_mydata1, y) #D = 0.048552, p-value = 0.9283
Fn = ecdf(log_mydata1)
Fm = ecdf(y)
plot(Fn,
main="Cumulative Distribution Functions",
xlab="data",
ylab="Cumulative Frequency",
pch=NA, lwd= 2,
col = "red")
lines(Fm, pch=NA, lty=1, lwd= 2)
rug(log_mydata1)
rug(y, side = 3, col = "red")
legend("topleft",
legend=c("original", "simulated"),
lty=c(1,1),
col=c("black", "red"))
grid()
[![enter image description here][3]][3]
The `p-value` of Kolmogorov-Smirnov test is `0.9283` and CDFs are close each to other.
But according to the documentaion of the [JohnsonDistribution package][4] I should use the **z**JohnsonDistribution function instead of the **y**JohnsonDistribution function.
Also I confused with `sim` series. In my case, `sim` is the normal distributed variable, but it is should [be uniformly distributed on the unit interval \[0, 1\]][1].
**Questions.** Am I correct in my steps? How to inverse correctly the Johnson normalized variable to a marginal variable? Should I use the `qnorm()` function?
**Edit 2.**
I have changed the library to do the direct/back Johnson transformation:
library(Johnson)
set.seed(1)
n <- 252
log_mydata1 <- runif(n, min=-0.06, max= 0.05)
# Applying SU Johnson transformation
jt_mydata1<-RE.Johnson(log_mydata1)
z1 <- jt_mydata1$transformed; shapiro.test(z1) # W = 0.99495, p-value = 0.5748
mean(z1); sd(z1) # -0.02276707, 1.002103
# fpar1 <- ArimaModelFit(z1)
# Coefficients:
# ar1 ar2 ar3 ma1 ma2 intercept
# 0.5520 0.4231 -0.1234 -0.5468 -0.4532 -0.0302
# s.e. 0.3582 0.3492 0.0666 0.3572 0.3570 0.0075
#ARMA Simulation
sim1 <- arima.sim(list(order = c(3,0,2),
ar = c(0.5520, 0.4231, -0.1234),
ma = c(-0.5468, -0.4532)), n = n)
mean(sim1);sd(sim1) # -0.0350542, 1.041598
shapiro.test(sim1) # W = 0.99285, p-value = 0.2675
# convert an normalized variable back to a marginal variable
gamma1 <- jt_mydata1$f.gamma
lambda1 <- jt_mydata1$f.lambda
epsilon1 <- jt_mydata1$f.epsilon
eta1 <- jt_mydata1$f.eta
inv_jt_mydata1 <- lambda1 * sinh((sim1 - gamma1)/eta1) + epsilon1
mean(inv_jt_mydata1);sd(inv_jt_mydata1)
#[1] -0.09473195
#[1] 0.4666171
The p-value of Kolmogorov-Smirnov test is less 2.2e-16 and
ks.test(log_mydata1, inv_jt_mydata1)
# D = 0.53571, p-value < 2.2e-16
# alternative hypothesis: two-sided
CDFs are different:
[![enter image description here][5]][5]
[1]: https://en.wikipedia.org/wiki/Johnson's_SU-distribution
[2]: https://www.quantstart.com/articles/ARIMA-GARCH-Trading-Strategy-on-the-SP500-Stock-Market-Index-Using-R
[3]: https://i.sstatic.net/WfSm0.png
[4]: https://cran.r-project.org/web/packages/JohnsonDistribution/JohnsonDistribution.pdf
[5]: https://i.sstatic.net/e3CZa.jpg