Quantile Matching using the skewed t-distribution from Azzalini & Capitanio (2003)

I try to replicate some findings from a paper (page 11/12 of https://www.newyorkfed.org/medialibrary/media/research/staff_reports/sr914.pdf). In this paper they estimated some quantiles using quantile regressions. I did this step, so I have the 5%, 25%, 75% and 95% quantiles for multiple periods and variables. Then, they used these estimated quantiles to fit a conditional probability distribution for every period and variable. The probability distribution used is the four-parameter skew t-distribution of Azzalini & Capitanio (2003). I found the R package, which generates this kind of distribution (https://cran.r-project.org/web/packages/sn/sn.pdf). So now I don’t know how to proceed. My code that I wrote is the following:

I have to use constraints since the omega and nu values can only be positive.

# the data which contains the quantiles(date: first column and 4 columns with quantiles (0.05, 0.25, 0.75, 0.95) for each period (192 rows) for one variable
imp_quantiles_rgdp<- read_excel("estimated_quantiles") 


density_range_rgdp = seq(from=-8,to= 8,by=0.1)

skewed_density_rgdp = matrix(ncol=length(density_range_rgdp), nrow=nrow(imp_quantiles_rgdp))

for (i in 1:nrow(imp_quantiles_rgdp)){
  imp_q <- as.numeric(imp_quantiles_rgdp[i,2:5]) 
        # the estimated quantiles generated in the first step using the quantile regression framework 

  q_fc <- function(x)sum((imp_q-qst(c(0.05,0.25,0.75,0.95),xi=x[1],omega=x[2],alpha=x[3],nu=x[4]))^2)
        # the idea is to minimize the squared difference between the skew t-implied quantiles and the estimated quantiles 

  osol <- optim(c(0,1,0,1),fn=q_fc, lower=c(-Inf, 0, -Inf, 0.0000001),  method="L-BFGS-B")
  skewed_para <- osol$par
  den <- dst(density_range_rgdp, xi=skewed_para[1], omega=skewed_para[2], alpha=skewed_para[3], nu=skewed_para[4])
  skewed_density_rgdp[i,] <- den

}

The problem I have now:

• It does not work for some variables. It runs forever, usually it gets stuck at the same row depending on the initial values. When I use different initial values, it sometime works, but then the results are totally wrong.

Depending on the initial values, I get errors like:

Error in qsn(p, xi, omega, alpha) :
failed convergence, try with solver="RFB"

Error in if (na == 1 & nz > 3 & all(alpha * z > -5) & (tau == 0L)) "T.Owen" else "biv.nt.prob" :
missing value where TRUE/FALSE needed
Called from: psn(x, dp = dp0, …)

So my guess is that there is an easier and more efficient way to solve my problem. Moreover, I guess that I often estimate local min instead of a global. I am thinking about using a different optimization method or maybe there is even a package, which matches quantiles. An additional concern was that my estimated quantiles are not always monotonically increasing, so I probably have to use an uncrossing procedure. Like the author of the paper mention on page 11 in the footnote. But this does not seem to be the issue, when the method gets stuck, since the quantiles in this row fulfill this condition. Furthermore, it works sometimes regardless of the non-increasing quantiles. I guess the non-increasing quantiles would explain the wrong/ weird distributions, but not why it get stuck.

Hopefully, someone can help me. Thanks.