Is it possible to use a normal copula after the Johnson transformation? I have a sample, log-returns of a finance time series. The sample size is n=252. I'd like to estimate a distribution in order to modeling dependence with a copula.
The binned histogram of frequencies, cumulative frequency counts, and qq-plot are below:

I have applied the normality hypothesis test (Anderson-Darling test) and the null hypothesis was rejected at level 0.05 (p-value is 0.003905).
I have applied the Kolmogorov-Smirnov test and the null hypothesis can't rejected at level 0.05 (p-value is 0.21769). But the null hypothesis should be rejected based on the Shapiro-Wilk test (p-value is 0.211e-4), the Chen-Shapiro test (5% critical value is 0.00356).
Original data is here CFC.csv
mydata = read.csv("CFC.csv")  # read csv file 

log_returns <- diff(log(mydata$Close), lag=1)
n <- length(log_returns)
mydata$log_ret = with(mydata, c(NA, diff(log(mydata$Close))))
#range(log_returns)
par(mfrow=c(1,3))
hist(log_returns, freq = TRUE, label=TRUE, breaks = 14, xlim=c(-0.08,0.08))
curve(dnorm(x, mean=mean(log_returns), 
               sd=sd(log_returns)), add=TRUE, col="darkblue", lwd=2)

Fn = ecdf(log_returns) 
plot(Fn, 
   main="CDF of log_returns", 
   xlab="log_returns", 
   ylab="Cumulative Frequency") 
qqnorm(log_returns)
qqline(log_returns)

shapiro.test(log_returns)
# Shapiro-Wilk normality test

# data:  log_returns
# W = 0.97509, p-value = 0.000211

# required 
library(nortest)
ad.test(log_returns)
# Anderson-Darling normality test

# data:  log_returns
# A = 1.1991, p-value = 0.003905

Problem. Data is not normal. Add I don't know how to specify a distribution to modeling dependence with a copula.
My attempt is: I have applied the Johnson's transformation to normalize the data.
# require
library(Johnson)
#Applying Johnson transformation 
log_returns_JT<-RE.Johnson(log_returns)

#$p
#[1] 0.2304269

#$f.gamma
#[1] 0.1002716

#$f.lambda
#[1] 0.03428872

#$f.epsilon
#[1] 0.00456466

#$f.eta
#[1] 2.109682

hist(log_returns_JT$transformed, freq = TRUE, label=TRUE)

curve(dnorm(x, mean=mean(log_returns), 
               sd=sd(log_returns)), add=TRUE, col="darkblue", lwd=2)

curve(dnorm(x, mean=mean(log_returns_JT$transformed), 
               sd=sd(log_returns_JT$transformed)), add=TRUE, col="red", lwd=2)

Fn = ecdf(log_returns_JT$transformed) 
plot(Fn, 
   main="CDF of log_returns_JT", 
   xlab="log_returns_JT", 
   ylab="Cumulative Frequency") 
qqnorm(log_returns_JT$transformed)
qqline(log_returns_JT$transformed)

shapiro.test(log_returns_JT$transformed)
#Shapiro-Wilk normality test

#data:  log_returns_JT$transformed
#W = 0.99465, p-value = 0.5217

ad.test(log_returns_JT$transformed)
#Anderson-Darling normality test

#data:  log_returns_JT$transformed
#A = 0.48099, p-value = 0.2304

(adt.p<-RE.ADT(log_returns_JT$transformed)$p)
# [1] 0.2304269

Here are the binned histogram of frequencies, cumulative frequency counts, and qq-plot after the Johnson's transformation:

I have applied the normality hypothesis test (Anderson-Darling test) and the null hypothesis can't reject at level 0.05 (p-value is 0.2304269).
Edit. I'm reading the papers recommended by @user25459 but, unfortunately, I don't find the answer for my question yet.
Questions. Is it possible to use a normal copula to modeling dependence and set the marginal distribution as the normal distribution? 
After the modeling I'm going to apply the inverse transformation.
 A: There is now a proof that derives the complete set of asset and liability classes, though all are mixture distributions.   You can find it at https://ssrn.com/abstract=2828744
The reason for the mixture is bankruptcy risk, merger risk, liquidity costs, and a stochastic budget constraint. Some debt instruments can be modeled using the normal or lognormal, but no equity-type asset can have a first moment.  Some assets, such as antiques, have the ratio of two Gumbel distributions as its distribution so there is no analytic representation. 
Generally, there is no distribution,  although for the broad population of US equities, its remarkably close to a truncated Cauchy distribution.  Multiply it by a very wide survival function and you'll get an extraordinarily good fit.   You can think of a trade as a survivor of the budget constraint because,  of course, not all orders fill.
The fundamental argument of the paper is that returns are not data, instead prices are data and returns are a  transformation of prices.  In particular,  it is a future value divided by a present value with an appraisal error upon entry and exit, minus one.  As such, all returns are ratio distrubutions.   Zero coupons differ in that the payoff is essentially fixed, so part of the ratio is a constant.
For completeness,  a population study was done on US trades from 1925-2013 to test the most important claims.  It can be found at https://ssrn.com/abstract=2653151
edit  I saw after my post it was log returns.  It will be a potentially ugly transformation of the hyperbolic secant distribution.   Really we should quit using the log approximation,  the only reason the process started was the need to use slide rules.
edit 2 I crudely mapped your bins in Excel and it was close to the hyperbolic secant, though probably skewed.  You can do the Johnson transformation because you can map any univariate standard distribution to any other univariate standard distribution, but you need to be careful that you do not forget that properties that govern are the properties of the untransformed data.  
Consider the simplest possible case, the unrealistic case of Markowitz.  The distribution of returns in a universe with only two equal sized equity assets would be $$\frac{1}{2\pi}\frac{\gamma}{(\gamma^2+(r_1-\mu_1)^2+(r_2-\mu_2)^2)^{1.5}}.$$  You could map this distribution to any bivariate distribution, but look closely at its properties.
$\gamma$ is not a vector, it is a scalar.  If you add dimensions $\gamma$ will become $\gamma'$, but $\gamma'$ is still a scalar.  The likelihood has spherical uncertainty and the density has spherical errors.  You can map this to a normal distribution with a covariance matrix, but it will not be unique because you are mapping to a distribution with more parameters than exist in the raw data.
Let's imagine instead you map a Cauchy to a normal, with zero covariance on the off diagonals, that is $\sigma_{i,j}\equiv{0},\forall{i,j}$ and equal variances to assure circularity, that is $\sigma_{i,i}\equiv\sigma_{j,j},\forall{i,j}$.  In the normal distribution, asset 1 and asset 2 are independent, but in the raw data they are not independent.  You have preserved the shape, but not the dynamics.  The Cauchy distribution is an example of a distribution where none of the dimensions are correlated, but none are independent.
There is no theoretical reason to not map the raw data to a copula function, but the "how" to do that is something I cannot answer.  I have done very light reading on the topic as I have thought about it myself since the Cauchy distribution doesn't really appear except in a mixture, but I read enough to realize I would have to spend a lot of time on the topic.
My preference to solve the problem, if I ever get to it, would be to note that things such as bankruptcy and liquidity costs do covary with underlying factors and so I would use Bayesian model selection and tools such as marginalization to separate out the components.
A: 
Questions. Is it possible to use a normal copula to modeling dependence and set the marginal distribution as the normal distribution?
After the modeling I'm going to apply the inverse transformation.

As a basic principle, it's certainly possible to transform margins to approximate normality and then fit a multivariate Gaussian, which is in effect fitting a Gaussian copula to the dependence relationship.

I don't know how to specify a distribution to modeling dependence with a copula.

Perhaps I misunderstand your intent here, I'm not sure; the point of a copula is to fit the dependence structure independent of the marginals. Usually in copula work this is achieved by transforming the marginals to approximate uniformity (via one of several different approaches) then fitting some copula, but in the case of the Gaussian copula that's not necessarily the easiest approach (nor even the most common).
I see no particular issue with using Johnson transformations if they're adequate to achieve very near normality. It shouldn't be worse than some of the other approaches people use.

However, for asset returns usually the Gaussian copula is regarded as unsuitable because the tail dependence in the Gaussian copula goes to 0, while in real asset returns can show substantial tail dependence (in particular, if the market goes south, like say in the GFC, correlated assets will tend to see very extreme results at the same time -- the Gaussian copula can't capture that.
Sometimes the t-copula is used to model this, since it does have tail dependence, but it doesn't capture the asymmetry (the dependence is not equally strong in the upper and lower tail). It may be adequate for some purposes, however.
