How to convert a Johnson normalized variable back to a marginal variable

Question

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.

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 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()

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 I should use the zJohnsonDistribution function instead of the yJohnsonDistribution 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].

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:

Edit 3.

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
# SU
# inv_jt_mydata1 <- lambda1 * sinh((sim1 - gamma1)/eta1) +  epsilon1
# SB
inv_jt_mydata1 <- ((lambda1 + epsilon1)* exp((sim1 - gamma1)/eta1) +  epsilon1)/(1+exp((sim1 - gamma1)/eta1))

The p-value of Kolmogorov-Smirnov test is $0.97$

ks.test(log_mydata1, inv_jt_mydata1)
#    D = 0.043651, p-value = 0.97
# alternative hypothesis: two-sided

and CDFs are close each to other:

Answer 1

Edit after question edit 2

From the docs and the evidence, the RE.Johnson function seems to produce standard Normal variates:

beta_variates <- rbeta(3000,2,3)
jt_variates <- RE.Johnson(beta_variates)
hist(jt_variates$transformed, breaks = 40, probability = T, 
   main = "Normal density versus Johnson-transformed Beta variates")
points( jt_variates$transformed, dnorm(jt_variates$transformed))

The ARIMA simulator can mimic these well, so I suspect that the difference in your CDFs is caused by the back-transform. Indeed, skipping the simulation and back-transforming the original data does not work.

inv_jt_mydata1 <-  jt_variates$f.lambda * 
sinh((jt_variates - jt_variates$f.gamma)/jt_variates$f.eta) +
  jt_variates$f.epsilon
mean(inv_jt_mydata1);sd(inv_jt_mydata1)
qqplot(log_mydata1, inv_jt_mydata1, main = "Q-Q plot of original versus back-transformed data") 

Based on browsing the source code and the fact that jt_variates[[2]] == "SB", I suspect that RE.Johnson has applied the logit-like function seen here in 1.2, not the sinh in 1.3.

Original answer

Here's how I'm reading your procedure. You:

• Generate an iid uniform time-series
• Pretend it's Johnson-distributed
• Apply a transform that converts Johnson variates to standard normal variates, and does who-knows-what to your uniform variates
• Fit an arima model to your (IID) data
• Simulate from an arima with the fitted parameters using normal errors (You didn't change the rand.gen arg in arima.sim; it defaults to standard normal.)
• Apply a transform that converts standard normal variates to Johnson variates and does who-knows-what to your arima output.

You want to know why the final result looks like your original data. It's probably because the particular Johnson distribution you've chosen does a great job mimicking your original unifs.

jd_forward_transform = function(x) yJohnsonDistribution(x, iType, gamma, delta, lambda , xi)
hist(jd_transform(rnorm(50000)), breaks = 40)

As you note via the Shapiro test, the Johnson-to-normal transform has turned your almost-Johnson unifs into almost-normals. These are almost-mimicked by the exactly normal marginal distribution of the ARMA variates. In brief, who-knows-what is approximately correct, because statistical software is distressingly robust.