R/dist-stableFit.R

In fBasics: Rmetrics - Markets and Basic Statistics

# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Library General Public License for more details.
#
# You should have received a copy of the GNU Library General
# Public License along with this library; if not, write to the
# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston,
# MA  02111-1307  USA


################################################################################
# FUNCTION:            STABLE DISTRIBUTION:
#  stableFit            Fits parameters of a stable density
#  .phiStable            Creates contour table for McCulloch estimators
#  .PhiStable            Contour table created by function .phiStable()
#  .qStableFit           Estimates parameters by McCulloch's approach
#  .mleStableFit         Estimates stable parameters by MLE approach
#  .stablePlot           Plots results of stable parameter estimates
# LIBRARY:             DESCRPTION:
#  stabledist           Stable Duistribution
################################################################################


stableFit <-
    function(x, alpha = 1.75, beta = 0, gamma = 1, delta = 0,
    type = c("q", "mle"), doplot = TRUE, control = list(),
    trace = FALSE, title = NULL, description = NULL)
{
    # A function implemented by Diethelm Wuertz

    # Description
    #   Stable Parameter Estimation

    # FUNCTION:

    # Start Values: Use Quantile Method:
    ans = .qStableFit(x, doplot, title, description)

    # Continue with MLE Approach:
    if (type[1] == "mle") {
        Alpha = ans@fit$estimate[1]
        Beta  = ans@fit$estimate[2]
        Gamma = ans@fit$estimate[3]
        Delta = ans@fit$estimate[4]
        if (is.na(Alpha)) Alpha = alpha
        if (is.na(Beta)) Beta = beta
        if (is.na(Gamma)) Gamma = gamma
        if (is.na(Delta)) Delta = delta
        ans = .mleStableFit(x, Alpha, Beta, Gamma, Delta, doplot,
            control, trace, title, description)
    }

    # Return Value:
    ans
}


# ------------------------------------------------------------------------------


.phiStable <- 
    function()
{
    # A function implemented by Diethelm Wuertz

    # Description:
    #   Creates contour table for McCulloch estimators

    # Note:
    #   Stable Distribution - delta=1 and gamma=0 fixed!

    # FUNCTION:

    # Settings:
    alpha = c(seq(0.50, 1.95, by = 0.1), 1.95, 1.99)
    beta = c(-0.95, seq(-0.90, 0.90, by = 0.10), 0.95)
    m = length(alpha)
    n = length(beta)

    # phi1:
    phi1 = function(alpha, beta)
    {
        ( stabledist::qstable(0.95, alpha, beta) - 
          stabledist::qstable(0.05, alpha, beta) ) /
        ( stabledist::qstable(0.75, alpha, beta) - 
          stabledist::qstable(0.25, alpha, beta) )
    }

    # phi2:
    phi2 = function(alpha, beta)
    {
        ( ( stabledist::qstable(0.95, alpha, beta) - 
            stabledist::qstable(0.50, alpha, beta) ) -
          ( stabledist::qstable(0.50, alpha, beta) - 
            stabledist::qstable(0.05, alpha, beta) ) ) /
          ( stabledist::qstable(0.95, alpha, beta) - 
            stabledist::qstable(0.05, alpha, beta) )
    }

    # Phi:
    Phi1 = Phi2 = matrix(rep(0, n*m), ncol = n)
    for ( i in 1:m ) {
        for ( j in 1:n ) {
            Phi1[i,j] = phi1(alpha[i], beta[j])
            Phi2[i,j] = phi2(alpha[i], beta[j])
            print( c(alpha[i], beta[j], Phi1[i,j], Phi2[i,j]) )
        }
    }

    #Plot:
    contour(alpha, beta, Phi1, levels = c(2.5, 3, 5, 10, 20, 40),
        xlab = "alpha", ylab = "beta", labcex = 1.5, xlim = c(0.5, 2.0))
    contour(alpha, beta, Phi2, levels = c(-0.8, -0.6, -0.4, -0.2, 0,
        0.2, 0.4, 0.6, 0.8), col = "red",  labcex = 1.5, add = TRUE)

    # Result:
    .PhiStable <- list(Phi1 = Phi1, Phi2 = Phi2, alpha = alpha, beta = beta)

    # Dump:
    if (FALSE) dump(".PhiStable", "PhiStable.R")

    # Return Value:
    .PhiStable
}


# ------------------------------------------------------------------------------


.PhiStable <-
    structure(

    #  Contour table created by .phiStable()
    list(

    Phi1 = structure(c(28.1355600962322, 15.7196640771722,
    10.4340722145276, 7.72099712337154, 6.14340919629241, 5.14081963876442,
    4.45927409495436, 3.97057677836415, 3.60450193720703, 3.31957820409758,
    3.09091185492284, 2.9029090263133, 2.74659551412658, 2.61782756430233,
    2.51560631141019, 2.47408801877228, 2.44525363752631, 28.3483300097440,
    15.8035063527296, 10.4727953966709, 7.74141823634036, 6.15576314177463,
    5.14915506532249, 4.46531988542058, 3.97505189771183, 3.60736868847651,
    3.32104702700423, 3.09111342044363, 2.90219825204304, 2.74552084086904,
    2.61707392900167, 2.51531961206501, 2.47401179006119, 2.44525268403340,
    28.9048147202382, 16.0601313263480, 10.6117049719320, 7.8237753182648,
    6.20955682299195, 5.18679961805621, 4.49216000186678, 3.99373973242068,
    3.61933010719953, 3.3275594162867, 3.09331226057454, 2.90155994459038,
    2.74368539137288, 2.61555847928520, 2.51479444686423, 2.47387109904253,
    2.44525093764076, 29.8026200973137, 16.5596243974224, 10.9360864131028,
    8.04563520952583, 6.35889612112486, 5.28397230078563, 4.55431822343012,
    4.03266384946799, 3.64246979110989, 3.33973615437707, 3.09810405845187,
    2.90191073816964, 2.7422687239301, 2.61421528380001, 2.51433181630644,
    2.47374639426837, 2.44524854839106, 31.1711919144382, 17.3256645589998,
    11.4257682914949, 8.37916374365985, 6.59111369941092, 5.4435297531592,
    4.65827423282770, 4.09593909348039, 3.67827370972420, 3.35788028768518,
    3.10538249512478, 2.90322430476313, 2.7411383952518, 2.61307358265729,
    2.51393064512464, 2.47363793686619, 2.44524473760244, 33.1948934984822,
    18.4071500759968, 12.0794579487667, 8.80181848446003, 6.8744516193679,
    5.6351929385469, 4.78568154439060, 4.17585254205654, 3.72404890551037,
    3.38077811930023, 3.11457663716798, 2.90522469731234, 2.74043557937936,
    2.61209710262770, 2.51359601849030, 2.47354590407541, 2.44523688699559,
    35.9578028398622, 19.7936882741385, 12.8664734285429, 9.28349135864713,
    7.18278427917665, 5.83729171895546, 4.91779143967590, 4.25960094973301,
    3.77307337255403, 3.4057709031586, 3.12469898537975, 2.90767817230573,
    2.7399388028832, 2.61130578115398, 2.51330169831474, 2.47347041290067,
    2.4452243360266, 38.9647018181828, 21.2269841758554, 13.6587928142908,
    9.75737930074977, 7.47910500668215, 6.02679068365276, 5.03975173429821,
    4.33641171958119, 3.81844774991854, 3.42948489228359, 3.13453596844362,
    2.91006724722144, 2.73965159366498, 2.61068126597691, 2.51309087992836,
    2.47342375837911, 2.44524639012562, 41.7818369588849, 22.4602363370435,
    14.3128567516975, 10.1435099877839, 7.72005687380391, 6.18048140792379,
    5.13779081982581, 4.39769218775219, 3.85488105109873, 3.44880791704519,
    3.14262443869137, 2.91215029990006, 2.73947335432536, 2.61026149741526,
    2.51293721037603, 2.47333463795113, 2.44523871005695, 43.8597008728084,
    23.308181530155, 14.7432127411421, 10.3930748979735, 7.87518470289541,
    6.2796283727218, 5.20095216574708, 4.43724315622136, 3.87838777226609,
    3.46128790559551, 3.14792806986726, 2.91355836113947, 2.7394048824397,
    2.60999397178590, 2.51287974154683, 2.4733261698261, 2.44523291863013,
    44.6351177921226, 23.612190997173, 14.8937659530932, 10.4791077126402,
    7.9285095195091, 6.31375867971341, 5.2228691158165, 4.45085115209764,
    3.88645952709398, 3.46562279514434, 3.14978921717977, 2.91402779433739,
    2.73939000980242, 2.60993524495366, 2.51289337121820, 2.47330017392183,
    2.44524525291192, 43.8597008728085, 23.3081815301550, 14.7432127411420,
    10.3930748979735, 7.87518470289542, 6.2796283727218, 5.20095216574709,
    4.43724315622136, 3.87838777226609, 3.46128790559551, 3.14792806986726,
    2.91355836113947, 2.7394048824397, 2.60999397178590, 2.51287974154683,
    2.4733261698261, 2.44523291863013, 41.7818369588849, 22.4602363370436,
    14.3128567516975, 10.1435099877840, 7.72005687380392, 6.18048140792379,
    5.13779081982581, 4.39769218775219, 3.85488105109874, 3.44880791704519,
    3.14262443869137, 2.91215029990006, 2.73947335432536, 2.61026149741526,
    2.51293721037603, 2.47333463795113, 2.44523871005695, 38.9647018181828,
    21.2269841758554, 13.6587928142909, 9.75737930074976, 7.47910500668216,
    6.02679068365277, 5.0397517342982, 4.33641171958119, 3.81844774991854,
    3.42948489228359, 3.13453596844362, 2.91006724722144, 2.73965159366498,
    2.61068126597691, 2.51309087992836, 2.47342375837911, 2.44524639012562,
    35.9578028398623, 19.7936882741384, 12.8664734285429, 9.28349135864712,
    7.18278427917667, 5.83729171895546, 4.91779143967590, 4.25960094973301,
    3.77307337255403, 3.4057709031586, 3.12469898537975, 2.90767817230572,
    2.7399388028832, 2.61130578115398, 2.51330169831474, 2.47347041290067,
    2.4452243360266, 33.1948934984822, 18.4071500759967, 12.0794579487667,
    8.80181848446, 6.8744516193679, 5.6351929385469, 4.78568154439059,
    4.17585254205653, 3.72404890551037, 3.38077811930023, 3.11457663716797,
    2.90522469731234, 2.74043557937936, 2.61209710262770, 2.51359601849030,
    2.47354590407541, 2.44523688699559, 31.1711919144381, 17.3256645589999,
    11.4257682914948, 8.37916374365983, 6.59111369941092, 5.44352975315921,
    4.65827423282769, 4.09593909348039, 3.67827370972420, 3.35788028768518,
    3.10538249512478, 2.90322430476312, 2.7411383952518, 2.61307358265729,
    2.51393064512464, 2.47363793686619, 2.44524473760245, 29.8026200973136,
    16.5596243974222, 10.9360864131028, 8.04563520952581, 6.35889612112487,
    5.28397230078562, 4.55431822343012, 4.03266384946799, 3.64246979110989,
    3.33973615437707, 3.09810405845187, 2.90191073816964, 2.7422687239301,
    2.61421528380001, 2.51433181630644, 2.47374639426837, 2.44524854839106,
    28.904814720238, 16.0601313263481, 10.6117049719320, 7.82377531826481,
    6.20955682299195, 5.18679961805620, 4.49216000186678, 3.99373973242068,
    3.61933010719952, 3.3275594162867, 3.09331226057455, 2.90155994459037,
    2.74368539137288, 2.61555847928520, 2.51479444686423, 2.47387109904253,
    2.44525093764076, 28.3483300097438, 15.8035063527296, 10.4727953966709,
    7.74141823634038, 6.15576314177465, 5.1491550653225, 4.46531988542057,
    3.97505189771183, 3.60736868847651, 3.32104702700424, 3.09111342044363,
    2.90219825204304, 2.74552084086904, 2.61707392900167, 2.51531961206501,
    2.47401179006119, 2.44525268403340, 28.1355600962322, 15.7196640771723,
    10.4340722145275, 7.72099712337151, 6.1434091962924, 5.14081963876442,
    4.45927409495435, 3.97057677836415, 3.60450193720703, 3.31957820409757,
    3.09091185492284, 2.9029090263133, 2.74659551412658, 2.61782756430233,
    2.51560631141019, 2.47408801877228, 2.44525363752631),
    .Dim = as.integer(c(17, 21))),

    Phi2 = structure(c(-0.984831521754346, -0.96178750287388,
    -0.925815277018118, -0.877909786944587, -0.820184465666658,
    -0.755031918545081, -0.684592921700777, -0.610570901457792,
    -0.534175805219629, -0.456236240234657, -0.377306837254648,
    -0.297867095171864, -0.218666685919607, -0.141143940005887,
    -0.0674987199592516, -0.0328701597273596, -0.00642990194045526,
    -0.984833720280198, -0.96124832221758, -0.92402252050786,
    -0.87419335313835, -0.81415246030398, -0.74662653505109,
    -0.674063099101783, -0.598388680797427, -0.520979416823626,
    -0.442739494396008, -0.364272193075272, -0.286088382456793,
    -0.208957354222478, -0.1342789733294, -0.0640144100029875,
    -0.0311490271069921, -0.006091726394748, -0.979153629182453,
    -0.95106168789293, -0.909673941239617, -0.85602448791855,
    -0.791828629762056, -0.720369142925194, -0.644748722301942,
    -0.567214164881932, -0.489279988980024, -0.411880943403832,
    -0.335605401995838, -0.260987626602553, -0.188786833912990,
    -0.120269390643364, -0.0570145517746281, -0.0277025095388298,
    -0.0054153433806392, -0.956099620177616, -0.91510100889139,
    -0.863494586547947, -0.804340589970157, -0.739811913292119,
    -0.670992067490536, -0.598508014750715, -0.523938483324115,
    -0.449119649196787, -0.375320171482799, -0.303334211640077,
    -0.233819259537015, -0.167632931803862,
    -0.105970645511563, -0.0499759998221647, -0.0242509183116938,
    -0.00473892214395457, -0.91091521430624, -0.854269872890404,
    -0.792380659464175, -0.728324327412516, -0.663641470500818,
    -0.598766118303959, -0.533390833007158, -0.466958108675312,
    -0.399718402638334, -0.332810169846175, -0.267433026630897,
    -0.204649764654678, -0.145582752525226, -0.0913961065988934,
    -0.042902827489912, -0.0207949382390745, -0.004062467796571,
    -0.838142578217344, -0.767730556629547, -0.699387193072487,
    -0.634346095207953, -0.572568000770375, -0.513428073110925,
    -0.456016815864829, -0.399237447140644, -0.342164056355244,
    -0.284798223846443, -0.228151639956042, -0.173628819918976,
    -0.122674652213850, -0.0765638152109311, -0.0358016086121756,
    -0.0173352283908653, -0.00338423605633909, -0.732960648615699,
    -0.655211639338905, -0.586717060962992, -0.525770503418058,
    -0.470570307907465, -0.419573011635011, -0.371434023729476,
    -0.324886123523844, -0.278750785874156, -0.232309518692667,
    -0.185933156371725, -0.140978784025337, -0.0990680854456507,
    -0.0615175707971035, -0.0286670029088599, -0.0138724303903631,
    -0.00270395845154647, -0.592691443906747, -0.517917221384721,
    -0.456902886405761, -0.405327421204478, -0.360251537992231,
    -0.319685291504111, -0.282167838420900, -0.246496132463733,
    -0.211596523966543, -0.176602222829403, -0.141397972727901,
    -0.106937786396409, -0.0748572933895385, -0.0462908166605395,
    -0.0215205174484398, -0.0104071747947459, -0.00203767706388774,
    -0.418561278354954, -0.359043286074154, -0.313110693536944,
    -0.275656631235694, -0.243720459128403, -0.215493360906447,
    -0.18976206222188, -0.165585707487220, -0.142157842427747,
    -0.118764989711438, -0.0951531051692878, -0.0718822809574559,
    -0.0501702371114745, -0.0309413187181704, -0.0143569748225324,
    -0.00693928545290556, -0.00135914851834483, -0.217058958479026,
    -0.183957984146554, -0.159253973812649, -0.139529858864502,
    -0.122956050946611, -0.108468261383110,
    -0.095363933833862, -0.0831489044941912, -0.0713807547092049,
    -0.0596749091239935, -0.0478406346829651, -0.0361189030171852,
    -0.0251638752545685, -0.0154951657479659, -0.00718251696063777,
    -0.00347157285038577, -0.000682830247122471, -1.30195010613456e-15,
    -7.12101008020212e-16, -9.8472751183435e-16, -1.13492047149954e-15,
    -1.21096413626332e-15, -1.05505102779748e-15, -1.11782289130157e-15,
    -8.13224377427617e-16, -5.85149540688739e-16, 4.61239028385592e-16,
    1.45510576814449e-16, -4.733934001018e-16, -3.36773719381392e-16,
    0, 9.23542032962673e-17, -1.87944839722067e-16, -2.8550539469612e-16,
    0.217058958479025, 0.183957984146551, 0.159253973812647,
    0.139529858864502, 0.12295605094661, 0.108468261383109,
    0.0953639338338615, 0.0831489044941909, 0.0713807547092042,
    0.0596749091239931, 0.0478406346829641, 0.036118903017185,
    0.0251638752545679, 0.0154951657479665, 0.00718251696063749,
    0.00347157285038558, 0.00068283024712209, 0.418561278354956,
    0.359043286074156, 0.313110693536942, 0.275656631235694,
    0.243720459128404, 0.215493360906444, 0.189762062221879,
    0.165585707487219, 0.142157842427746, 0.118764989711438,
    0.0951531051692875, 0.0718822809574556, 0.0501702371114741,
    0.0309413187181703, 0.0143569748225322, 0.00693928545290556,
    0.00135914851834492, 0.59269144390675, 0.517917221384721,
    0.456902886405764, 0.405327421204478, 0.360251537992232,
    0.319685291504111, 0.2821678384209, 0.246496132463732,
    0.211596523966543, 0.176602222829402, 0.141397972727901,
    0.106937786396409, 0.074857293389538, 0.0462908166605391,
    0.0215205174484394, 0.0104071747947454, 0.00203767706388774,
    0.7329606486157, 0.655211639338903, 0.58671706096299,
    0.525770503418055, 0.470570307907467, 0.419573011635009,
    0.371434023729475, 0.324886123523844, 0.278750785874155,
    0.232309518692667, 0.185933156371724, 0.140978784025336,
    0.0990680854456504, 0.0615175707971032, 0.0286670029088598,
    0.0138724303903627, 0.00270395845154656, 0.838142578217343,
    0.767730556629544, 0.699387193072486, 0.634346095207951,
    0.572568000770374, 0.513428073110924, 0.456016815864827,
    0.399237447140644, 0.342164056355244, 0.284798223846443,
    0.228151639956042, 0.173628819918976, 0.122674652213850,
    0.0765638152109311, 0.0358016086121755, 0.0173352283908648,
    0.00338423605633909, 0.910915214306239, 0.854269872890404,
    0.792380659464175, 0.728324327412515, 0.663641470500817,
    0.598766118303959, 0.533390833007157, 0.466958108675312,
    0.399718402638334, 0.332810169846175, 0.267433026630897,
    0.204649764654678, 0.145582752525226, 0.0913961065988931,
    0.0429028274899118, 0.020794938239074, 0.00406246779657053,
    0.956099620177616, 0.915101008891389, 0.863494586547946,
    0.804340589970156, 0.739811913292119, 0.670992067490536,
    0.598508014750714, 0.523938483324114, 0.449119649196786,
    0.375320171482798, 0.303334211640076, 0.233819259537015,
    0.167632931803862, 0.105970645511562, 0.0499759998221642,
    0.0242509183116936, 0.00473892214395428, 0.979153629182453,
    0.95106168789293, 0.909673941239617, 0.85602448791855,
    0.791828629762055, 0.720369142925193, 0.644748722301942,
    0.567214164881931, 0.489279988980024, 0.411880943403833,
    0.335605401995838, 0.260987626602552, 0.188786833912991,
    0.120269390643364, 0.0570145517746277, 0.0277025095388299,
    0.00541534338063891, 0.984833720280198, 0.96124832221758,
    0.92402252050786, 0.87419335313835, 0.81415246030398,
    0.746626535051091, 0.674063099101783, 0.598388680797427,
    0.520979416823626, 0.442739494396008, 0.364272193075271,
    0.286088382456793, 0.208957354222478, 0.1342789733294,
    0.0640144100029874, 0.0311490271069919, 0.00609172639474771,
    0.984831521754346, 0.96178750287388, 0.925815277018118,
    0.877909786944586, 0.820184465666658, 0.75503191854508,
    0.684592921700776, 0.610570901457792, 0.534175805219629,
    0.456236240234657, 0.377306837254648, 0.297867095171864,
    0.218666685919607, 0.141143940005886, 0.0674987199592513,
    0.0328701597273591, 0.00642990194045459),
    .Dim = as.integer(c(17, 21))),

    alpha = c(0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.1, 1.2, 1.3, 1.4,
        1.5, 1.6, 1.7, 1.8, 1.9, 1.95, 1.99),

    beta = c(-0.95, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1,
        0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95)),

    .Names = c("Phi1", "Phi2", "alpha", "beta")
)


# ------------------------------------------------------------------------------


.qStableFit <- 
    function(x, doplot = TRUE, title = NULL, description = NULL)
{
    # A function implemented by Diethelm Wuertz

    # Description:
    #   Estimates stable parameters by McCulloch's approach

    # Note:
    #   This implementation assumes delta=1 and gamma=0

    # FUNCTION:

    # Settings:
    CALL = match.call()

    # Load Contour Table:
    # data(PhiStable)
    Phi1 = .PhiStable$Phi1
    Phi2 = .PhiStable$Phi2
    alpha = .PhiStable$alpha
    beta = .PhiStable$beta

    # Estimate phi:
    r = sort(x)
    N = length(r)
    q95 = r[round(0.95*N)]
    q75 = r[round(0.75*N)]
    q50 = r[round(0.50*N)]
    q25 = r[round(0.25*N)]
    q05 = r[round(0.05*N)]
    phi1 = max( (q95-q05) / (q75-q25), 2.4389 )
    phi2 = ((q95-q50)-(q50-q05)) / (q95-q05)
    # print(c(min(Phi1), phi1, max(Phi1)))
    # print(c(min(Phi2), phi2, max(Phi2)))

    # Plot:
    if (doplot) {
        contour(alpha, beta, Phi1, levels = c(2.5, 3, 5, 10, 20, 40),
            xlab = "alpha", ylab = "beta", xlim = c(0.5, 2.0))
        contour(alpha, beta, Phi2, levels = c(-0.8, -0.6, -0.4, -0.2,
            0.2, 0.4, 0.6, 0.8), col = "red", add = TRUE)
        lines(c(0.5, 2), c(0, 0), col = "red")
        contour(alpha, beta, Phi1, levels = phi1, add = TRUE, lty = 3,
            col = "blue")
        contour(alpha, beta, Phi2, levels = phi2, add = TRUE, lty = 3,
            col = "blue")
        title(main = "Stable Quantiles")
    }

    # Extract Estimate from Contours, if possible:
    u = contourLines(alpha, beta, Phi1, levels = phi1)
    Len = length(u)
    if( Len > 0) {
        u = u[[1]][-1]
        v = contourLines(alpha, beta, Phi2, levels = phi2)
        # print("v")
        # print(v)
        v = v[[1]][-1]
        xout = seq(min(v$y), max(v$y), length = 200)
        z = approx(v$y, v$x, xout = xout)$y - approx(u$y, u$x, xout = xout)$y
        index = which.min(abs(z))
        V = round(xout[index], 3)
        U = round(approx(v$y, v$x, xout = xout[index])$y, 3)
        if (doplot) points(U, V, pch = 19, cex = 1)
    } else {
        U = V = NA
    }

    # Add Title and Description:
    if (is.null(title)) title = "Stable Parameter Estimation"
    if (is.null(description)) description = ""

    if (is.na(U) | is.na(V)) {
        GAM = NA
    } else {
        phi3 = stabledist::qstable(0.75, U, V) - 
               stabledist::qstable(0.25, U, V)
        GAM = (q75-q25) / phi3
    }

    if (is.na(U) | is.na(V)) {
        DELTA = NA
    } else {
        phi4 = -stabledist::qstable(0.50, U, V) + V*tan(pi*U/2)
        DELTA = phi4*GAM - V*GAM*tan(pi*U/2) + q50
    }

    # Fit:
    fit = list(estimate = c(alpha = U, beta = V, gamma = GAM, delta = DELTA))

    # Return Value:
    new("fDISTFIT",
        call = as.call(CALL),
        model = "Stable Distribution",
        data = as.data.frame(x),
        fit = fit,
        title = as.character(title),
        description = description )
}


# ------------------------------------------------------------------------------


.mleStableFit <- 
    function(x, alpha = 1.75, beta = 0, gamma = 1, delta = 0, doplot = TRUE,
    control = list(), trace = FALSE, title = NULL, description = NULL)
{
    # A function implemented by Diethelm Wuertz

    # Description:
    #   Estimates stable parameters by MLE approach

    # Note:
    #   This implementation assumes delta=1 and gamma=0

    # FUNCTION:

    # Transform:
    x.orig = x
    x = as.vector(x)

    # Settings:
    CALL = match.call()

    # Log-likelihood Function:
    obj = function(x, y = x, trace = FALSE) {
        f = -mean(log(stabledist::dstable(y, 
            alpha = x[1], beta = x[2], gamma = x[3], delta = x[4])))
        # Print Iteration Path:
        if (trace) {
            cat("\n Objective Function Value:  ", -f)
            cat("\n Stable Estimate:           ", x)
            cat("\n")
        }
        f
    }
    
    # Minimization:
    eps = 1e-4
    r <- nlminb(
        objective = obj, 
        start = c(alpha, beta, gamma, delta), 
        lower = c( eps, -1+eps, 0+eps, -Inf),
        upper = c(2-eps, 1-eps,  Inf,  Inf),
        y = x, control = control, trace = trace) 
    alpha = r$par[1]
    beta = r$par[2]
    gamma = r$par[3]
    delta = r$par[4]
    
    # Optional Plot:
    if (doplot) .stablePlot(x, alpha, beta, gamma, delta)

    # Add Title and Description:
    if (is.null(title)) title = "Stable Parameter Estimation"
    if (is.null(description)) description = ""

    # Fit:
    fit = list(estimate = c(alpha, beta, gamma, delta), 
        minimum = -r$objective, code = r$convergence, gradient = r$gradient)

    # Return Value:
    new("fDISTFIT",
        call = as.call(CALL),
        model = "Stable Distribution",
        data = as.data.frame(x.orig),
        fit = fit,
        title = as.character(title),
        description = description )
}


# ------------------------------------------------------------------------------


.stablePlot <-
    function(x, alpha, beta, gamma, delta)
{
    # A function implemented by Diethelm Wuertz

    # Description:
    #   Estimates stable parameters by MLE approach
    
    # FUNCTION:
    
    # Plot:
    span.min <- stabledist::qstable(0.01, alpha, beta, gamma, delta)
    span.max <- stabledist::qstable(0.99, alpha, beta, gamma, delta)
    span = seq(span.min, span.max, length = 100)
    par(err = -1)
    z = density(x, n = 100)
    x = z$x[z$y > 0]
    y = z$y[z$y > 0]
    y.points = stabledist::dstable(span, alpha, beta, gamma, delta)
    ylim = log(c(min(y.points), max(y.points)))
    plot(x, log(y), xlim = c(span[1], span[length(span)]),
        ylim = ylim, type = "p", xlab = "x", ylab = "log f(x)")
    title("Stable Distribution: Parameter Estimation")
    lines(x = span, y = log(y.points), col = "steelblue")
    grid()
    
    # Return Value:
    invisible()
}


################################################################################