R/methods-PME.R

In pmoments: Partial Moments

#################################################################################
##
##   R package pmoments by Alexios Ghalanos Copyright (C) 2008
##   This file is part of the R package pmoments.
##
##   The R package pmoments is free software: you can redistribute it and/or modify
##   it under the terms of the GNU General Public License as published by
##   the Free Software Foundation, either version 3 of the License, or
##   (at your option) any later version.
##
##   The R package pmoments 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 General Public License for more details.
##
#################################################################################

pmExpectation<-function(data, f, method="sample", fit=NULL, threshold=0, moment=0, type="lower", standardize=FALSE, 
		analytical = "quadrature", n.boot=1000, subdivisions=1000, tailThresh= 10, data.names=NULL)
{
	UseMethod("pmExpectation")
}


.pmExpectation<-function(data, f, method="sample", fit=NULL, threshold=0, moment=0, type="lower", standardize=FALSE, 
		analytical = "quadrature", n.boot=1000, subdivisions=1000, tailThresh= 10, data.names=NULL)
{
	# check data is available or distributions
	if(!missing(data)){
		nc=dim(as.matrix(data))[2]
		if(nc==1) datatype="single" else datatype="multiple"
		if(length(method)==1) method=rep(method,nc)
		if(is.null(data.names)){
			if(is.null(colnames(as.matrix(data)))) data.names=paste("X",1:nc,sep="") else data.names=colnames(as.matrix(data))
		}
	} else{
		nc=length(f)
		if(nc==1) datatype="single" else datatype="multiple"
		data=as.matrix(NA)
		if(is.null(data.names)) data.names=paste("X",1:nc,sep="")
		if(length(data.names)!=nc) data.names=paste("X",1:nc,sep="")
		if(length(method)==1) method=rep(method,nc)
	}
	if(missing(f)) f=NA
	# case1: univariate data
	# case2: multivariate data
	ans<-switch(datatype,
			single=.pmExpectationI(data,f,method,fit,threshold,moment,type,standardize,analytical,n.boot,subdivisions, tailThresh, data.names),
			multiple=.pmExpectationII(data,f,method,fit,threshold,moment,type,standardize,analytical,n.boot,subdivisions, tailThresh, data.names))
	return(ans)
}

setMethod("pmExpectation",signature(data="matrix",f="vector"),.pmExpectation)
setMethod("pmExpectation",signature(data="matrix",f="missing"),.pmExpectation)
setMethod("pmExpectation",signature(data="missing",f="vector"),.pmExpectation)
setMethod("pmExpectation",signature(data="vector",f="vector"),.pmExpectation)
setMethod("pmExpectation",signature(data="vector",f="missing"),.pmExpectation)


.pmExpectationII<-function(data, f, method="sample", fit, threshold, moment=0, type="lower", standardize=FALSE, 
		analytical = "quadrature", n.boot=1000, subdivisions=1000, tailThresh= 10, data.names)
{
	# check method dimension
	n<-length(data.names)
	if(length(method)==1) method = rep(method,n)
	pmsol<-vector(mode="list",length=n)
	xfit<-vector(mode="list",length=n)
	distribution<-vector(mode="character",length=n)
	# check for mixed (analytical&sample) methods
	if(any(method=="analytical")){
		z<-which(method=="analytical")
		zn<-length(z)
		xf<-vector(mode="character",length=n)
		if(length(f)>zn) xf[z]<-f[z] else xf[z]<-f
		if(length(analytical)!=1 && length(analytical)!=zn) stop("analytical option length does not match corresponding method length")
		if(length(f)>zn){
			for(i in 1:length(z)){
				xfit[[z[i]]]<-fit[[z[i]]]
				distribution[z]<-f[z]
			}
		} else {
			for(i in 1:length(z)){
			xfit[[z[i]]]<-fit[[i]]
			}
			distribution[z]<-f
		}
		if(length(analytical)==1) xanalytical=rep(analytical,zn) else xanalytical=analytical
		zxanalytical<-vector(mode="character",length=n)
		zxanalytical[z]<-xanalytical
	}
	# checks on the sample/boot/jack methods and mixed methods
	if(missing(data) & any(method=="sample" | method=="boot" | method=="jack")) stop("cannot calculate sample partial moments without a data object!")
	if(!missing(data) & any(method=="sample" | method=="boot" | method=="jack")){
		s<-which(method=="sample" | method=="boot" | method=="jack")
		sdata<-as.matrix(data)
		distribution[s]<-"empirical"
	}
	# n is either the dimension of the dataset or the length of the distribution vector
	if(length(threshold)==1) xthreshold=rep(threshold,n) else xthreshold=threshold
	# technically, mixing moments is bad as the object cannot be then used in portfolio
	# optimization...use pmoments instead for univariate results.
	if(length(moment)!=1){
		if(length(moment)!=n) stop("Moment vector length not equal to data width. Either supply the correct length or single moment for all data.") else xmoment=moment
	} else{
	xmoment=rep(moment,n)
	}
	if(length(type)==1) xtype=rep(type,n) else xtype=type
	# again mixing this is bad.
	if(length(standardize)==1) xstandardize=rep(standardize,n) else xstandardize=standardize
	if(length(subdivisions)==1) xsubdivisions=rep(subdivisions,n) else xsubdivisions=subdivisions	
	if(length(tailThresh)==1) xtailThresh=rep(tailThresh,n) else xtailThresh=tailThresh
	# dispatch to method function
	for(i in 1:n){
		pmsol[[i]]<-switch(method[i],
				analytical=.pmAnalytical(f=xf[i], threshold=xthreshold[i], fit=xfit[[i]], subdivisions=xsubdivisions[i], moment=xmoment[i], 
						tailThresh=xtailThresh[i], type=xtype[i], standardize=xstandardize[i], analytical=zxanalytical[i]),
				sample=.pmSample(sdata[,i], xthreshold[i], xmoment[i], xtype[i], xstandardize[i]),
				boot = .pmBootStrap(sdata[,i], n.boot=n.boot, xthreshold[i], xmoment[i], xtype[i], xstandardize[i], returnType="II"),
				jack = .pmJackKnife(sdata[,i], xthreshold[i], xmoment[i], xtype[i], xstandardize[i],returnType="II"))
	}
	sol<-new("PME",datanames=data.names,
			standardize = xstandardize,
			tail=type,
			moment = xmoment,
			method = method,
			threshold = xthreshold,
			distribution = distribution,
			pm = pmsol,
			data = data,
			datatype="multiple")
	return(sol)
}

.pmExpectationI<-function(data, f, method="sample", fit, threshold, moment=c(0,1,2,3,4), 
		type="lower", standardize=FALSE, analytical = "quadrature", n.boot=1000, subdivisions=1000, tailThresh= 10, data.names)
{	
	# must allow multiple distributions.
	# check for mixed (analytical&sample) methods
	nm=length(method)
	n = length(moment)
	if(!missing(f)) nf=length(f) else nf=0
	distribution=vector(mode="character",length=(nm-1) + nf)
	if(any(method=="analytical")){
		zf=which(method=="analytical")	
		distribution<-rep(NA,(nm-1)+nf)
		distribution[zf:(zf+nf-1)]<-f
		distribution[which(is.na(distribution))]<-"empirical"
		if(is.na(f[1])) stop("cannot have analytical method with missing distribution name(f)")
		if(is.null(fit)) stop("cannot have analytical without distribution fit object")
		#if(!is.na(f[1]) & length(f)!=length(fit)) stop("length of f not equal to length of fit")
	}
	if(missing(data) & (any(method=="sample") | any(method=="boot") | any(method=="jack"))) stop("cannot calculate sample partial moments without a data object!")
	if(!missing(data) & (any(method=="sample") | any(method=="boot") | any(method=="jack"))){
		if(any(method!="analytical")){
		distribution[which(method=="sample" | method=="boot" | method=="jack")]<-"empirical"
	}}
	pmsol = vector(mode="list",length=nm)
	names(pmsol)<-method
	# dispatch to method function
	for(j in 1:nm){
		msol=vector(mode="list",length=n)
		names(msol)<-paste("moment",round(moment,3),sep="")
		for(i in 1:n)
		{
			msol[[i]]<-switch(method[j],
				analytical=.pmAnalyticalMultiple(f=f, threshold=threshold, fit=fit, subdivisions=subdivisions, moment=moment[i], 
						tailThresh=tailThresh, type=type, standardize=standardize, analytical=analytical),
				sample=.pmSample(as.numeric(data), threshold, moment[i], type, standardize),
				boot = .pmBootStrap(x=as.numeric(data), n.boot=n.boot, threshold, moment[i], type, standardize, returnType="II"),
				jack = .pmJackKnife(x=as.numeric(data), threshold, moment[i], type, standardize, returnType="II"))}
	pmsol[[j]]=msol
	}
	# must have something in the data slot (no matrixOrNULL class)
	sol<-new("PME",datanames=data.names,
			standardize = standardize,
			tail=type,
			moment = moment,
			method = method,
			threshold = threshold,
			distribution = distribution,
			pm = pmsol,
			data = as.matrix(data),
			datatype="single")	
	return(sol)
}

#----------------------------------------------------------------------------------
# Partial moments estimation functions
# Analytical (quadrature & trapezium)
# Sample Based
# Boostrap and JackKnife
# Consider Robust?

.pmAnalyticalMultiple<-function(f, threshold, fit, subdivisions=1000, moment=0, tailThresh= 10, 
		type="lower", standardize=FALSE, analytical="quadrature")
{
	n=length(f)
	pasol<-vector(mode="numeric",length=n)
	names(pasol)<-f
	for(i in 1:n){
		pasol[i]<-.pmAnalytical(f=f[i], threshold=threshold, fit=fit[[i]], subdivisions=subdivisions, moment=moment, tailThresh= tailThresh, 
				type=type, standardize=standardize, analytical=analytical)
	}
	return(pasol)
}

.pmAnalytical<-function(f, threshold, fit, subdivisions=1000, moment=0, tailThresh= 10, 
		type="lower", standardize=FALSE, analytical="quadrature")
{
	ans<-switch(analytical,
			quadrature=.pmIntegrate(f, threshold, fit, subdivisions, moment, tailThresh, type, standardize),
			trapezium=.pmNumIntegrate(f,threshold,fit, subdivisions, moment, tailThresh, type, standardize))
	return(ans)
}
# Partial Moments Quadrature Integration method
.pmIntegrate <- function(f, threshold, fit, subdivisions=1000, moment=0, tailThresh= 10, 
		type="lower", standardize=FALSE)
{
	#reflection: lower and upper
	if(type == "lower")
	{
		sgn = -1
		tailThresh = sgn*abs(tailThresh)
		lower = tailThresh
		upper = threshold
	}
	else
	{
		sgn = 1
		tailThresh = sgn*abs(tailThresh)
		lower = threshold
		upper = tailThresh
	}
	if(standardize) m = moment else m = 1
	#f <- match.fun(f)
	ff <- eval(parse(text=paste("function(x) d",f,"(x,fit) * (sgn*(-threshold + x))^moment",sep="")))
	# some pdf functions are problemtic on the etxremes (-Inf or large subdivisions)
	ans = try(expr=integrate(ff, lower, upper,subdivisions=subdivisions,stop.on.error = FALSE), silent = TRUE)
	if(inherits(ans, "try-error")){
		ans<-integrate(ff, lower, upper,subdivisions=5,stop.on.error = FALSE)}
	else{ans<-integrate(ff, lower, upper,subdivisions=subdivisions,stop.on.error = FALSE)}
	if(moment==0) ans=ans$value else ans=ans$value^(1/m)
	return(ans)
}

# Partial Moments Trapezium Integration method
.pmNumIntegrate <- function(f,threshold,fit, subdivisions=1000, moment=0, tailThresh= 10, 
		type="lower", standardize=FALSE)
{
	if(type == "lower")
	{
		sgn = -1
		tailThresh = sgn*abs(tailThresh)		
		lower = tailThresh
		upper = threshold
	}
	else
	{
		sgn = 1
		tailThresh = sgn*abs(tailThresh)
		lower = threshold
		upper = tailThresh
	}
	if(standardize) m = moment else m = 1
	#f = match.fun(f)
	ff <- eval(parse(text=paste("function(x) d",f,"(x,fit)",sep="")))
	Ri = seq(lower, upper, length.out=subdivisions)
	dR = mean(diff(Ri))
	pRi = ff(Ri)
	d1 = .sumC(pRi)/2
	d2 = d1*dR
	sR = .sumC((sgn*(-threshold+Ri))^moment)/2
	if(moment==0) ans=sum(d2*sR) else ans=sum(d2*sR)^(1/m)
	return(ans)
}

# Partial Moments Sample Method
.pmSample<-function(x, threshold, moment, type="lower", standardize=FALSE)
{
	pm = switch(type,
			lower = .lpmSample(x, threshold, moment),
			upper = .upmSample(x, threshold, moment))
	if(standardize) m = moment else m = 1
	if(moment==0) ans=pm else ans=pm^(1/m)
	names(pm)<-"sample"
	return(ans)
}

# Lower Sample Method
.lpmSample<-function(x,threshold,moment)
{
	if(moment==0)
	{
		ans = sum(as.integer(x<=threshold))/length(x)
	}
	else
	{
		ans = (mean(.max.linear(threshold-x)^moment))
	}
	return(ans)
}

# Upper Sample Method
.upmSample<-function(x,threshold,moment)
{
	if(moment==0)
	{
		ans = mean(as.integer(x>threshold))
	}
	else
	{
		ans = (mean(.max.linear(x-threshold)^moment))
	}
	return(ans)
}

# Boostrap and JackKnife
.pmBootStrap<-function(x, n.boot=1000, threshold, moment, type="lower", standardize=FALSE, 
		returnType="I")
{
	bt = bootstrap(x=x, nboot=n.boot,theta=.pmSample,threshold=threshold, moment=moment, 
			type=type, standardize=standardize, func=NULL)
	if(returnType=="I") ans = mean(bt$theta) else ans = bt
	return(ans)
}

.pmJackKnife<-function(x, threshold, moment, type="lower", standardize=FALSE, 
		returnType="I")
{
	jk = jackknife(x=x, theta=.pmSample,threshold=threshold, moment=moment, 
			type=type, standardize=standardize)
	if(returnType=="I") ans = mean(jk$jack.values) else ans = jk
	return(ans)
}

.pmValues<-function(x, threshold, moment, type="lower", standardize=FALSE)
{
	pm = switch(type,
			lower = .lpmValues(x, threshold, moment),
			upper = .upmValues(x, threshold, moment))
	if(standardize) m = moment else m = 1
	if(moment==0) ans=pm else ans=pm^(1/m)
	names(pm)<-"sample"
	return(ans)
}

# Lower Sample Method
.lpmValues<-function(x,threshold,moment)
{
	if(moment==0)
	{
		ans = cumsum((as.integer(x<=threshold))/length(x))
	}
	else
	{
		ans = (cumsum(.max.linear(threshold-x)^moment))/length(x)
	}
	return(ans)
}

# Upper Sample Method
.upmValues<-function(x,threshold,moment)
{
	if(moment==0)
	{
		ans = cumsum((as.integer(x>threshold))/length(x))
	}
	else
	{
		ans = (cumsum(.max.linear(x-threshold)^moment))/length(x)
	}
	return(ans)
}