R/fitStable.R
In jansila/mscPack:
Defines functions
quantFit
mleFit
print.myStable
summary.myStable
phiFun
dyn.load("src/stable.so")
# Core of functions implemented by Diethelm Wuertz, minor adjustments by Jan Sila
quantFit <-function(x, title = NULL, doplot=FALSE,description = NULL){
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)
# 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
}
if (is.na(U) | is.na(V)) {
GAM = NA
} else {
phi3 = quantStable(0.75, U, V) -
quantStable(0.25, U, V)
GAM = (q75-q25) / phi3
}
if (is.na(U) | is.na(V)) {
DELTA = NA
} else {
phi4 = -quantStable(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))
class(fit)<-"stableFit"
# Return Value:
ans<-list("call" = call,
"model" = "Stable Distribution",
"fit" = fit)
class(ans)<-"myStable"
return(ans)
}
# ------------------------------------------------------------------------------
mleFit <-function(input, init=NULL, trace = FALSE, title = NULL, description = NULL){
call<-match.call()
if(is.null(init)){init=quantFit(input)$fit$estimate
x=as.vector(init)
alpha=x[1]
beta=x[2]
gamma=x[3]
delta=x[4]
}else{x=init}
if(any(is.na(x))){
x<-c(1.8,0,0.7,0)
names(x)<-c("alpha","beta","gamma","delta")
alpha=x[1]
beta=x[2]
gamma=x[3]
delta=x[4]
}
obj = function(input, x , trace = FALSE) {
# cat("\nparameters are: ", c(x[1],x[2],x[3],x[4],"\n"))
f = -mean(log(denStable(input,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")
}
return(f)
}
# Minimization:
eps = 1e-4
r <- nlminb(
objective = obj,
start = c(x[1], x[2], x[3], x[4]),
lower = c( 1.5+2*eps, -0.95+eps, 0+eps, -Inf),
upper = c(2-eps, 0.95-eps, 1-eps, Inf),
input = input, trace = trace)
alpha = r$par[1]
beta = r$par[2]
gamma = r$par[3]
delta = r$par[4]
# Fit:
fit = list(estimate = c(alpha, beta, gamma, delta),
minimum = -r$objective, AIC=8-2*log(abs(r$objective)),code = r$convergence, gradient = r$gradient)
class(fit)<-"stableFit"
# Return Value:
ans<-list("call"=call, "fit"=fit)
class(ans)<-"myStable"
return(ans)
}
print.myStable<-function(object){
cat("Call: ");print(object$call)
cat("Estimate is: \n");print(object$fit$estimate)
cat("AIC: ");print(object$fit$AIC)
}
summary.myStable<-function(object){
cat("Call: ",object$call)
cat("Estimate is: ", object$fit$estimate)
cat("AIC: ",object$fit$AIC)
}
phiFun <-
function()
{
# Description:
# Creates contour table for McCulloch estimators
# 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)
{
( qstable(0.95, alpha, beta) -
qstable(0.05, alpha, beta) ) /
( qstable(0.75, alpha, beta) -
qstable(0.25, alpha, beta) )
}
# phi2:
phi2 = function(alpha, beta)
{
( ( qstable(0.95, alpha, beta) -
qstable(0.50, alpha, beta) ) -
( qstable(0.50, alpha, beta) -
qstable(0.05, alpha, beta) ) ) /
( qstable(0.95, alpha, beta) -
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")
)
# ------------------------------------------------------------------------------
################################################################################
jansila/mscPack documentation built on May 18, 2019, 2:38 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions quantFit mleFit print.myStable summary.myStable phiFun
dyn.load("src/stable.so")
# Core of functions implemented by Diethelm Wuertz, minor adjustments by Jan Sila
quantFit <-function(x, title = NULL, doplot=FALSE,description = NULL){
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)
# 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
}
if (is.na(U) | is.na(V)) {
GAM = NA
} else {
phi3 = quantStable(0.75, U, V) -
quantStable(0.25, U, V)
GAM = (q75-q25) / phi3
}
if (is.na(U) | is.na(V)) {
DELTA = NA
} else {
phi4 = -quantStable(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))
class(fit)<-"stableFit"
# Return Value:
ans<-list("call" = call,
"model" = "Stable Distribution",
"fit" = fit)
class(ans)<-"myStable"
return(ans)
}
# ------------------------------------------------------------------------------
mleFit <-function(input, init=NULL, trace = FALSE, title = NULL, description = NULL){
call<-match.call()
if(is.null(init)){init=quantFit(input)$fit$estimate
x=as.vector(init)
alpha=x[1]
beta=x[2]
gamma=x[3]
delta=x[4]
}else{x=init}
if(any(is.na(x))){
x<-c(1.8,0,0.7,0)
names(x)<-c("alpha","beta","gamma","delta")
alpha=x[1]
beta=x[2]
gamma=x[3]
delta=x[4]
}
obj = function(input, x , trace = FALSE) {
# cat("\nparameters are: ", c(x[1],x[2],x[3],x[4],"\n"))
f = -mean(log(denStable(input,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")
}
return(f)
}
# Minimization:
eps = 1e-4
r <- nlminb(
objective = obj,
start = c(x[1], x[2], x[3], x[4]),
lower = c( 1.5+2*eps, -0.95+eps, 0+eps, -Inf),
upper = c(2-eps, 0.95-eps, 1-eps, Inf),
input = input, trace = trace)
alpha = r$par[1]
beta = r$par[2]
gamma = r$par[3]
delta = r$par[4]
# Fit:
fit = list(estimate = c(alpha, beta, gamma, delta),
minimum = -r$objective, AIC=8-2*log(abs(r$objective)),code = r$convergence, gradient = r$gradient)
class(fit)<-"stableFit"
# Return Value:
ans<-list("call"=call, "fit"=fit)
class(ans)<-"myStable"
return(ans)
}
print.myStable<-function(object){
cat("Call: ");print(object$call)
cat("Estimate is: \n");print(object$fit$estimate)
cat("AIC: ");print(object$fit$AIC)
}
summary.myStable<-function(object){
cat("Call: ",object$call)
cat("Estimate is: ", object$fit$estimate)
cat("AIC: ",object$fit$AIC)
}
phiFun <-
function()
{
# Description:
# Creates contour table for McCulloch estimators
# 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)
{
( qstable(0.95, alpha, beta) -
qstable(0.05, alpha, beta) ) /
( qstable(0.75, alpha, beta) -
qstable(0.25, alpha, beta) )
}
# phi2:
phi2 = function(alpha, beta)
{
( ( qstable(0.95, alpha, beta) -
qstable(0.50, alpha, beta) ) -
( qstable(0.50, alpha, beta) -
qstable(0.05, alpha, beta) ) ) /
( qstable(0.95, alpha, beta) -
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")
)
# ------------------------------------------------------------------------------
################################################################################
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.