R/methods-plot.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.
##
#################################################################################

#----------------------------------------------------------------------------------
# Partial Moments Optimization plots
# Weights Plot
weightsPlot = function(object, col = NULL, legend = TRUE, title=NULL)
{
	UseMethod("weightsPlot")
}

.weightsPlot=function(object, col = NULL, legend = TRUE, title=NULL)
{   
	if (is.null(col)) col = rainbow(ncol(object@weights))
	weights = object@weights
	pos.weights = +0.5 * (abs(weights) + weights)
	neg.weights = -0.5 * (abs(weights) - weights)
	ymax = max(rowSums(pos.weights))
	ymin = min(rowSums(neg.weights))
	range = ymax - ymin
	ymax = ymax + 0.005 * range
	ymin = ymin - 0.005 * range
	dim = dim(weights)
	range = dim[1]
	xmin = 0
	xmax = range + 0.2 * range
	if(!legend){
		barplot(t(pos.weights), space = 0, ylab = "",
				ylim = c(ymin, ymax), col = col, border = "grey")
	} else {
		legendtext = colnames(object@weights)
		if(is.null(legendtext)){
			for(i in 1:dim[2]){legendtext[i] = paste("Asset", i, sep = " ")}
		}
		barplot(t(pos.weights), space = 0, ylab = "", xlim = c(xmin, xmax),
				ylim = c(ymin, ymax), col = col, border = "grey")
		legend("topright", legend = legendtext, bty = "n", cex = 0.8,
				fill = col)
	}
	barplot(t(neg.weights), space = 0, add = TRUE, col = col, border = "grey") 
	targetRisk = object@riskMeasure
	targetReturn = object@rewardMeasure 
	targetRaR = targetReturn/targetRisk
	nSigma = length(targetRisk)
	nLabels = 6
	M = c(0, ( 1:(nSigma %/% nLabels) ) ) *nLabels + 1
	nPrecision = 3
	axis(1, at = M, labels = signif(targetRisk[M], nPrecision),cex.axis=0.8)
	axis(3, at = M, labels = signif(targetReturn[M], nPrecision),cex.axis=0.8)
	
	mtext("Target Risk", side = 1, line = 2, cex = 0.9)
	mtext("Target Return", side = 3, line = 2, cex = 0.9)
	mtext("Weights (sum)", side = 2, line = 2, cex = 0.9)
	
	lines(x = c(0, nSigma), c(1, 1), col = "grey", lty = 3) 
	lines(x = c(0, nSigma), c(0, 0), col = "grey", lty = 3)   
	
	minIndex = which.min(targetRisk)
	minRisk = signif(min(targetRisk), 3)    
	# Add vertical Line at max risk/return:
	tR = which(targetRaR==max(targetRaR))
	tRR = signif(targetRisk[tR], 3)
	abline(v = tR, col = "steelblue", lty = 2, lwd = 2)
	
	mtext(paste(
					object@solverType, "|", object@solver, "|", "minRisk =", minRisk),
			side = 4, adj = 0, col = "grey", cex = 0.7)
	mtext(paste(
					object@solverType, "|", object@solver, "|", "optimalRisk =", tRR),
			side = 4, adj = 0, padj=1.2 , col = "grey", cex = 0.7)
	
	if(is.null(title)) mtext("Weights", adj = 0, line = 2.5, font = 2, cex = 0.8) else mtext(title, adj = 0, line = 2.5, font = 2, cex = 0.8)
	box()
	invisible()
}

setMethod("weightsPlot", signature(object="PMSOLVER"), .weightsPlot)

weightsPie = function(object, legend = TRUE, title=NULL)
{
	UseMethod("weightsPie")
}

.weightsPie<-function(object, legend = TRUE, title=NULL)
{     
	if(length(object@riskMeasure)!=1) stop("function not available for frontier")
	weights<-(object@weights)
	wsort<-sort.int(weights,index.return = TRUE)
	weights=weights[1,wsort$ix]
	labelsx = names(weights)
	labels = names(weights)
	Radius = 1
	if(length(weights)>15){
		z=which(abs(weights)<0.01)
		labels[z]=""
		Radius = 0.85
	}
	negw<-which(weights<0)
	posw<-which(weights>=0)
	negc<-.colorgradient(n = length(negw), low.col = 0.0, high.col=0.1, saturation = 1)
	posc<-.colorgradient(n = length(posw), low.col = 0.6, high.col=0.65, saturation = 1)
	mycol=vector(mode="character",length=length(weights))
	mycol=c(negc,posc)
	par(no.readonly = TRUE)
    pie(abs(weights), labels = labels, col = mycol, radius = Radius,cex=0.6)
    box()
    # Add Title:
	if(is.null(title)){
    title(main = paste("Partial Moment Optimization Weights [moment: ", object@moment,"]",sep=""))
	} else{
		title(main=title)
	}
    # Add Info:
    mtext(paste("Solver Type :", object@solverType, " | reward measure :", round(object@rewardMeasure,5)," | risk measure :",round(object@riskMeasure,5),sep=""), 
        side = 4, adj = 0, col = "grey", cex = 0.7)
	mtext(paste("Blue Tones: Positive Weights | Red Tones: Negative Weights",sep=""), 
		side = 1, adj=0, col = "grey", cex = 0.7)
	mtext(paste("Assets with abs(weight)<1% do not display labels",sep=""), 
		side = 1, adj=0, padj=1.5, col = "grey", cex = 0.7)
	mtext(paste("pmoments",sep=""), 
		side = 2, adj=0,padj=-1, col = "grey", cex = 0.7)
	if(legend){
        legendWeights = as.character(round(100*weights, digits = 1))
        legendWeights = paste(legendWeights, "%",sep=" ")
		legend(1.1,1, legend = labelsx, bty = "n", cex = 0.8, 
				fill = mycol)
        legend(1.3,1, legend = legendWeights,bty = "n", cex = 0.8)
    }
    # Return Value:
    invisible()
}

setMethod("weightsPie", signature(object="PMSOLVER"), .weightsPie)

# Frontier Plots & Lines
.plot.pmSolver<-function(x,...)
{
	if(length(x@riskMeasure)==1) stop("x is not a solver type IV (frontier) object")
	risk<-x@riskMeasure
	reward<-x@rewardMeasure
	riskFree=x@riskFree
	tangent = x@tangent
	xlabel<-paste("Lower Partial Moment [m = ",x@moment,"]",sep="")
	ylabel<-"Expected Return"
	minrisk = which(risk==min(risk))
	plot(risk[1:(minrisk-1)],reward[1:(minrisk-1)],type="p",col="black",ylab=ylabel,xlab=xlabel,xlim=c(0.5*min(risk),1.1*max(risk)),
			ylim=c(0,1.05*max(reward)),lwd=2,...)
	points(risk[minrisk:length(risk)],reward[minrisk:length(reward)],col="grey",lwd=2)
	points(risk[tangent],reward[tangent],col="steelblue",lwd=3)
	points(risk[minrisk],reward[minrisk],col="violet",lwd=3)
	title("LPM Frontier Plot",col="black",lty=2)
	leg.txt=c("Frontier(>minRisk)", "Frontier(<minrisk)",paste("Tangency(riskFree=",round(riskFree,4),")",sep=""), "minRisk")
	legend(x="bottomright", yjust=0,legend=leg.txt,lwd=3, lty=c(1,1,1,3), col=c('black','grey','steelblue','violet'), cex=0.8)
	mtext(paste("pmoments package"), side = 4, adj = 0, padj=0, col = "gray", cex = 0.7)
	grid()
	invisible()
}

.lines.pmSolver<-function(x,...)
{
	if(length(x@riskMeasure)==1) stop("x is not a solver type IV (frontier) object")
	risk<-x@riskMeasure
	reward<-x@rewardMeasure
	lines(risk,reward,...)
	invisible()
}

.points.pmSolver<-function(x,...)
{
	if(length(x@riskMeasure)==1) stop("x is not a solver type III (frontier) object")
	risk<-x@riskMeasure
	reward<-x@rewardMeasure
	points(risk,reward,...)
	invisible()
}

setMethod("plot", signature(x="PMSOLVER",y="missing"), .plot.pmSolver)
setMethod("lines", signature(x="PMSOLVER"), .lines.pmSolver)
setMethod("points", signature(x="PMSOLVER"), .points.pmSolver)

#----------------------------------------------------------------------------------
# Partial Moments Utility Plots
.plot.pmu<-function(x,...)
{
	dx<-sort(x@data)
	dy<-x@utility
	plot(dx,dy,type="l",ylab="Utility",xlab="X",...)
	leg.text<-paste("LM=",x@lmoment,", UM=",x@umoment,", LC=",x@lrc,", UC=",x@urc,sep="")
	#legend("topleft",legend=leg.text,text.col=,fill=1,cex=0.7,inset=c(0.03,0.06))
	title("Upper-to-Lower Partial Moment Utility Function")
	grid()
	invisible()
}


.lines.pmu<-function(x,...)
{
	op <- options()
	options(warn = -1)
	dx<-sort(x@data)
	dy<-x@utility
	lines(dx,dy,...)
	options(op)
	invisible()
}

setMethod("plot", signature(x="PMU",y="missing"), .plot.pmu)
setMethod("lines", signature(x="PMU"), .lines.pmu)

#----------------------------------------------------------------------------------
# Partial Moments Expectation Plots
# Sub-methods for single/multiple data type and estimation method
.plot.pme<-function(x,...)
{
	if(x@datatype=="single"){
	mm=x@method
	for(i in 1:length(mm)){
	sol<-switch(mm[i],
			boot=.sbootplot(x,...),
			jack=.sjackplot(x,...),
			sample=.ssampleplot(x,...),
			analytical=.sanalyticalplot(x,...))
	
	}} else{
	munq=rev(unique(x@method))
	for(i in 1:length(munq)){
	sol<-switch(munq[i],
			boot=.mbootplot(x,...),
			jack=.mjackplot(x,...),
			sample=.msampleplot(x,...),
			analytical=.manalyticalplot(x,...))
	}}
	invisible()
}

setMethod("plot", signature(x="PME",y="missing"), .plot.pme)

#----------------------------------------------------------------------------------
# Single DataType Plots
.sbootplot<-function(x,...)
{
	n = length(x@moment)
	cl=rainbow(n=n)
	if(n>16){
		scr=floor(n/16)
		z=n-16*floor(n/16)
		start=dev.next(which = dev.cur())
		tally=0
		for(j in 1:scr){
			dev.new(start+j)
			par(mfrow=c(4,4))
			for(i in 1:16){
				tmp=x@pm$boot[[i+tally]]$thetastar
				hist(tmp,col=cl[i+tally],ylab="frequency",xlab=paste("moment = ",round(x@moment[i+tally],2),sep=""),main="")
				mtext(paste("pmoments  : bootstrap"), side = 4, adj = 0, padj=0, col = "gray", cex = 0.4)
			}
			title(paste("Boostrapped Partial Moment Histograms [",x@datanames,"]","\n(page...",j,")",sep=""), outer=TRUE, line=-2, cex=0.75)
			tally=tally+16
		}
		if(z!=0){
			dev.new(dev.next(which = dev.cur())+1)
			par(mfrow=c(4,4))
			for(i in 1:z){
				tmp=x@pm$boot[[i+tally]]$thetastar
				hist(tmp,col=cl[i+tally],ylab="frequency",xlab=paste("moment = ",round(x@moment[i+tally],2),sep=""),main="")
				mtext(paste("pmoments  : bootstrap"), side = 4, adj = 0, padj=0, col = "gray", cex = 0.4)
			}
			title(paste("Boostrapped Partial Moment Histograms [",x@datanames,"]","\n(page...",j,")",sep=""), outer=TRUE, line=-2, cex=0.75)
		}
	} else{
		d=.divisortable(n)
		par(mfrow=c(d[1],d[2]))
		for(i in 1:n){
			tmp=x@pm$boot[[i]]$thetastar
			hist(tmp,col=cl[i],ylab="frequency",xlab=paste("moment = ",round(x@moment[i],2),sep=""),main="")
			mtext(paste("pmoments  : bootstrap"), side = 4, adj = 0, padj=0, col = "gray", cex = 0.5)
			mtext(paste("mean: ",round(mean(tmp),3*(as.integer(x@moment[i]/2)))," | sd: ",round(sd(tmp),3*(as.integer(x@moment[i]/2))),sep=""),
			side = 4, adj = 0, padj=1.2, col = "gray", cex = 0.5)
			title(paste("Boostrapped Partial Moment Histograms [",x@datanames,"]",sep=""),outer=TRUE,cex=0.85, line=-2)
		}
	}
	invisible()
}

.sjackplot<-function(x,...)
{
	n = length(x@moment)
	cl=rainbow(n=n)
	if(n>16){
		scr=floor(n/16)
		z=n-16*floor(n/16)
		start=dev.next(which = dev.cur())
		tally=0
		for(j in 1:scr){
			dev.new(start+j)
			par(mfrow=c(4,4))
			for(i in 1:16){
				tmp=x@pm$jack[[i+tally]]$jack.values
				hist(tmp,col=cl[i+tally],ylab="frequency",xlab=paste("moment = ",round(x@moment[i+tally],2),sep=""),main="")
				mtext(paste("pmoments  : JackKnife"), side = 4, adj = 0, padj=0, col = "gray", cex = 0.4)
			}
			title(paste("JackKnife Partial Moment Histograms [",x@datanames,"]","\n(page...",j,")",sep=""), outer=TRUE, line=-2, cex=0.75)
			tally=tally+16
		}
		if(z!=0){
			dev.new(dev.next(which = dev.cur())+1)
			par(mfrow=c(4,4))
			for(i in 1:z){
				tmp=x@pm$jack[[i+tally]]$jack.values
				hist(tmp,col=cl[i+tally],ylab="frequency",xlab=paste("moment = ",round(x@moment[i+tally],2),sep=""),main="")
				mtext(paste("pmoments  : JackKnife"), side = 4, adj = 0, padj=0, col = "gray", cex = 0.4)
			}
			title(paste("JackKnife Partial Moment Histograms [",x@datanames,"]","\n(page...",j,")",sep=""), outer=TRUE, line=-2, cex=0.75)
		}
	} else{
		d=.divisortable(n)
		par(mfrow=c(d[1],d[2]))
		for(i in 1:n){
			tmp=x@pm$jack[[i]]$jack.values
			hist(tmp,col=cl[i],ylab="frequency",xlab=paste("moment = ",round(x@moment[i],2),sep=""),main="")
			mtext(paste("pmoments  : JackKnife"), side = 4, adj = 0, padj=0, col = "gray", cex = 0.6)
			mtext(paste("JackKnife s.e.(",round(x@pm$jack[[i]]$jack.se,5),")",sep=""), side = 4, adj = 0, padj=1.2, col = "gray", cex = 0.6)
			title(paste("JackKnife Partial Moment Histograms [",x@datanames,"]",sep=""),outer=TRUE,cex=0.85, line=-2)
		}
	}
	invisible()
}

.ssampleplot<-function(x,...)
{
	n = length(x@moment)
	cl=rainbow(n=n)
	if(n>16){
		scr=floor(n/16)
		z=n-16*floor(n/16)
		start=dev.next(which = dev.cur())
		tally=0
		tmpdata=x@data
		threshold=x@threshold
		tmpsdata=sort(tmpdata)
		for(j in 1:scr){
			dev.new(start+j)
			par(mfrow=c(4,4))
			for(i in 1:16){
				tmp=.pmValues(tmpsdata, threshold=threshold, moment=x@moment[i+tally], type=x@tail, standardize=x@standardize)
				if(x@tail=="lower"){
					pd=max(which(tmpsdata<=threshold))
					zt=1:pd
					plot(tmpsdata[zt],tmp[zt],col=cl[i+tally],ylab=paste("moment = ",round(x@moment[i+tally],2)), xlab="data",main="",cex.lab=0.9)
				} else{
					pd=min(which(tmpsdata>threshold))
					zt=pd:length(tmpsdata)
					plot(tmpsdata[zt],tmp[rev(zt)],col=cl[i+tally],ylab=paste("moment = ",round(x@moment[i+tally],2)), xlab="data",main="",cex.lab=0.9)
				}
				plot(tmpsdata[zt],tmp[zt],col=cl[i+tally],ylab=paste("moment = ",round(x@moment[i+tally],2)), xlab="data",main="",cex.lab=0.9)
				mtext(paste("pmoments  : sample"), side = 4, adj = 0, padj=0, col = "gray", cex = 0.4)
			}
			title(paste("Cumulative Sample Partial Moment-to-Threshold [",x@datanames,"]","\n Tail=",x@tail," (page...",j,")",sep=""), outer=TRUE, line=-2, cex=0.75)
			tally=tally+16
		}
		if(z!=0){
			dev.new(dev.next(which = dev.cur())+1)
			par(mfrow=c(4,4))
			for(i in 1:z){
				tmp=.pmValues(tmpsdata, threshold=threshold, moment=x@moment[i+tally], type=x@tail, standardize=x@standardize)
				if(x@tail=="lower"){
					pd=max(which(tmpsdata<=threshold))
					zt=1:pd
					plot(tmpsdata[zt],tmp[zt],col=cl[i+tally],ylab=paste("moment = ",round(x@moment[i+tally],2)), xlab="data",main="",cex.lab=0.9)
				} else{
					pd=min(which(tmpsdata>threshold))
					zt=pd:length(tmpsdata)
					plot(tmpsdata[zt],tmp[rev(zt)],col=cl[i+tally],ylab=paste("moment = ",round(x@moment[i+tally],2)), xlab="data",main="",cex.lab=0.9)
				}
				mtext(paste("pmoments  : sample"), side = 4, adj = 0, padj=0, col = "gray", cex = 0.4)
			}
			title(paste("Cumulative Sample Partial Moment-to-Threshold [",x@datanames,"]","\n Tail=",x@tail," (page...",j,")",sep=""), outer=TRUE, line=-2, cex=0.75)
		}
	} else{
		d=.divisortable(n)
		tmpdata=x@data
		threshold=x@threshold
		tmpsdata=sort(tmpdata)
		par(mfrow=c(d[1],d[2]))
		for(i in 1:n){
			tmp=.pmValues(tmpsdata, threshold=threshold, moment=x@moment[i], type=x@tail, standardize=x@standardize)
			if(x@tail=="lower"){
				pd=max(which(tmpsdata<=threshold))
			    zt=1:pd
				plot(tmpsdata[zt],tmp[zt],col=cl[i],ylab=paste("moment = ",round(x@moment[i],2)), xlab="data",main="",cex.lab=0.9)
			} else{
				pd=min(which(tmpsdata>threshold))
				zt=pd:length(tmpsdata)
				plot(tmpsdata[zt],tmp[rev(zt)],col=cl[i],ylab=paste("moment = ",round(x@moment[i],2)), xlab="data",main="",cex.lab=0.9)
			}
			mtext(paste("pmoments  : sample"), side = 4, adj = 0, padj=0, col = "gray", cex = 0.4)
			title(paste("Cumulative Sample Partial Moment-to-Threshold [",x@datanames,"]\n Tail=",x@tail,sep=""),outer=TRUE,cex=0.85, line=-2)
		}
	}
	invisible()
}

.sanalyticalplot<-function(x,...)
{
	# sample versus analytical
	# 2 plots (density and sample with partial moment shading on threshold
	# and barplot with overlayed moments (sample vs analytical)
	
	print("single dataset/analytical method: no plot method for this type of data")
}


#----------------------------------------------------------------------------------
.mbootplot<-function(x,...)
{
	# we need to distinguish the methods:
	zz=which(x@method=="boot")
	n = length(zz)
	datanames=x@datanames[zz]
	cl=rainbow(n=n)
	if(n>16){
		scr=floor(n/16)
		z=n-16*floor(n/16)
		#start=dev.new(dev.next(which = dev.cur())+1)
		tally=0
		for(j in 1:scr){
			dev.new(dev.next(which = dev.cur())+1)
			par(mfrow=c(4,4))
			for(i in 1:16){
				tmp=x@pm[[zz[i+tally]]]$thetastar
				hist(tmp,col=cl[zz[i+tally]],ylab="frequency",xlab=paste("moment = ",round(x@moment[zz[i+tally]],2),sep=""),main=datanames[zz[i+tally]])
				mtext(paste("pmoments  : bootstrap"), side = 4, adj = 0, padj=0, col = "gray", cex = 0.4)
			}
			title(paste("Boostrapped Partial Moment Histograms (page...",j,")",sep=""), outer=TRUE, line=-1, cex=0.75)
			tally=tally+16
		}
		if(z!=0){
			dev.new(dev.next(which = dev.cur())+1)
			par(mfrow=c(4,4))
			for(i in 1:z){
				tmp=x@pm[[zz[i+tally]]]$thetastar
				hist(tmp,col=cl[zz[i+tally]],ylab="frequency",xlab=paste("moment = ",round(x@moment[zz[i+tally]],2),sep=""),main=datanames[zz[i+tally]])
				mtext(paste("pmoments  : bootstrap"), side = 4, adj = 0, padj=0, col = "gray", cex = 0.4)
			}
			title(paste("Boostrapped Partial Moment Histograms (page...",scr+1,")",sep=""), outer=TRUE, line=-1, cex=0.75)
		}
	} else{
		d=.divisortable(n)
		dev.new(dev.next(which = dev.cur())+1)
		par(mfrow=c(d[1],d[2]))
		for(i in 1:n){
			tmp=x@pm[[zz[i]]]$thetastar
			hist(tmp,col=cl[i],ylab="frequency",xlab=paste("moment = ",round(x@moment[zz[i]],2),sep=""),main=datanames[i])
			mtext(paste("pmoments  : bootstrap"), side = 4, adj = 0, padj=0, col = "gray", cex = 0.5)
			mtext(paste("mean: ",round(mean(tmp),3*(as.integer(x@moment[zz[i]]/2)))," | sd: ",round(sd(tmp),3*(as.integer(x@moment[zz[i]]/2))),sep=""),
					side = 4, adj = 0, padj=1.2, col = "gray", cex = 0.5)
			title("Boostrapped Partial Moment Histograms", outer=TRUE,cex=0.85, line=-1)
		}
	}
}

.mjackplot<-function(x,...)
{
	# we need to distibuish the methods:
	zz=which(x@method=="jack")
	n = length(zz)
	datanames=x@datanames[zz]
	cl=rainbow(n=n)
	if(n>16){
		scr=floor(n/16)
		z=n-16*floor(n/16)
		#start=dev.next(which = dev.cur())
		tally=0
		for(j in 1:scr){
			dev.new(dev.next(which = dev.cur())+1)
			par(mfrow=c(4,4))
			for(i in 1:16){
				tmp=x@pm[[zz[i+tally]]]$jack.values
				hist(tmp,col=cl[zz[i+tally]],ylab="frequency",xlab=paste("moment = ",round(x@moment[zz[i+tally]],2),sep=""),main=datanames[zz[i+tally]])
				mtext(paste("pmoments  : JackKnife"), side = 4, adj = 0, padj=0, col = "gray", cex = 0.4)
			}
			title(paste("JackKnife Partial Moment Histograms (page...",j,")",sep=""), outer=TRUE, line=-1, cex=0.75)
			tally=tally+16
		}
		if(z!=0){
			dev.new(dev.next(which = dev.cur())+1)
			par(mfrow=c(4,4))
			for(i in 1:z){
				tmp=x@pm[[zz[i+tally]]]$jack.values
				hist(tmp,col=cl[zz[i+tally]],ylab="frequency",xlab=paste("moment = ",round(x@moment[zz[i+tally]],2),sep=""),main=datanames[zz[i+tally]])
				mtext(paste("pmoments  : JackKnife"), side = 4, adj = 0, padj=0, col = "gray", cex = 0.4)
			}
			title(paste("JackKnife Partial Moment Histograms (page...",scr+1,")",sep=""), outer=TRUE, line=-1, cex=0.75)
		}
	} else{
		d=.divisortable(n)
		dev.new(dev.next(which = dev.cur())+1)
		par(mfrow=c(d[1],d[2]))
		for(i in 1:n){
			tmp=x@pm[[zz[i]]]$jack.values
			hist(tmp,col=cl[zz[i]],ylab="frequency",xlab=paste("moment = ",round(x@moment[zz[i]],2),sep=""),main=datanames[zz[i+tally]])
			mtext(paste("pmoments  : JackKnife"), side = 4, adj = 0, padj=0, col = "gray", cex = 0.5)
			mtext(paste("mean: ",round(mean(tmp),3*(as.integer(x@moment[zz[i]]/2)))," | sd: ",round(sd(tmp),3*(as.integer(x@moment[zz[i]]/2))),sep=""),
					side = 4, adj = 0, padj=0.8, col = "gray", cex = 0.5)
			mtext(paste("JackKnife s.e.(",round(x@pm[[zz[i]]]$jack.se,5),")",sep=""), side = 4, adj = 0, padj=1.2, col = "gray", cex = 0.6)
			title("JackKnife Partial Moment Histograms", outer=TRUE,cex=0.85, line=-2)
		}
	}
}
#----------------------------------------------------------------------------------
# used in choosing the mfrow
.divisortable<-function(n)
{
	z=matrix(c(1,1,1,
			   2,2,1,
			   3,2,2,
			   4,2,2,
			   5,2,3,
			   6,2,3,
			   7,2,4,
			   8,2,4,
			   9,3,3,
			   10,3,4,
			   11,3,4,
			   12,3,4,
			   13,4,4,
			   14,4,4,
			   15,4,4,
			   16,4,4),ncol=3,byrow=TRUE)
	d=which(n==z[,1])
	return(z[d,2:3])
}
#----------------------------------------------------------------------------------
.msampleplot<-function(x,...)
{
	# we need to distinguish the methods:
	zz=which(x@method=="sample")
	n = length(zz)
	datanames=x@datanames[zz]
	cl=rainbow(n=n)
	if(n>16){
		scr=floor(n/16)
		z=n-16*floor(n/16)
		#start=dev.new(dev.next(which = dev.cur())+1)
		tally=0
		for(j in 1:scr){
			dev.new(dev.next(which = dev.cur())+1)
			par(mfrow=c(4,4))
			for(i in 1:16){
				tmpdata=x@data[,zz[i+tally]]
				threshold=x@threshold[zz[i+tally]]
				tmpsdata=sort(tmpdata)
				tmp=.pmValues(tmpsdata, threshold=threshold, moment=x@moment[zz[i+tally]], type=x@tail[1], standardize=x@standardize[zz[i+tally]])
				if(x@tail[1]=="lower"){
					pd=max(which(tmpsdata<=threshold))
					zt=1:pd
					plot(tmpsdata[zt],tmp[zt],col=cl[zz[i+tally]],ylab=paste("moment = ",round(x@moment[zz[i+tally]],2)), xlab="data",main=datanames[i+tally],cex.lab=0.9)
				} else{
					pd=min(which(tmpsdata>threshold))
					zt=pd:length(tmpsdata)
					plot(tmpsdata[zt],tmp[rev(zt)],col=cl[zz[i+tally]],ylab=paste("moment = ",round(x@moment[zz[i+tally]],2)), xlab="data",main=datanames[i+tally],cex.lab=0.9)
				}
				mtext(paste("pmoments  : bootstrap"), side = 4, adj = 0, padj=0, col = "gray", cex = 0.4)
			}
			title(paste("Cumulative Sample Partial Moment-to-Threshold (page...",j,")",sep=""), outer=TRUE, line=-1, cex=0.75)
			tally=tally+16
		}
		if(z!=0){
			dev.new(dev.next(which = dev.cur())+1)
			par(mfrow=c(4,4))
			for(i in 1:z){
				tmpdata=x@data[,zz[i+tally]]
				threshold=x@threshold[zz[i+tally]]
				tmpsdata=sort(tmpdata)
				tmp=.pmValues(tmpsdata, threshold=threshold, moment=x@moment[zz[i+tally]], type=x@tail[1], standardize=x@standardize[zz[i+tally]])
				if(x@tail[1]=="lower"){
					pd=max(which(tmpsdata<=threshold))
					zt=1:pd
					plot(tmpsdata[zt],tmp[zt],col=cl[zz[i+tally]],ylab=paste("moment = ",round(x@moment[zz[i+tally]],2)), xlab="data",main=datanames[i+tally],cex.lab=0.9)
				} else{
					pd=min(which(tmpsdata>threshold))
					zt=pd:length(tmpsdata)
					plot(tmpsdata[zt],tmp[rev(zt)],col=cl[zz[i+tally]],ylab=paste("moment = ",round(x@moment[zz[i+tally]],2)), xlab="data",main=datanames[i+tally],cex.lab=0.9)
				}
				mtext(paste("pmoments  : bootstrap"), side = 4, adj = 0, padj=0, col = "gray", cex = 0.4)
				}
			title(paste("Cumulative Sample Partial Moment-to-Threshold (page...",scr+1,")",sep=""), outer=TRUE, line=-1, cex=0.75)
		}
	} else{
		d=.divisortable(n)
		dev.new(dev.next(which = dev.cur())+1)
		par(mfrow=c(d[1],d[2]))
		for(i in 1:n){
			tmpdata=x@data[,i]
			threshold=x@threshold[i]
			tmpsdata=sort(tmpdata)
			tmp=.pmValues(tmpsdata, threshold=threshold, moment=x@moment[i], type=x@tail[1], standardize=x@standardize[i])
			if(x@tail[1]=="lower"){
				pd=max(which(tmpsdata<=threshold))
				zt=1:pd
				plot(tmpsdata[zt],tmp[zt],col=cl[i],ylab=paste("moment = ",round(x@moment[i],2)), xlab="data",main=datanames[i],cex.lab=0.9)
			} else{
				pd=min(which(tmpsdata>threshold))
				zt=pd:length(tmpsdata)
				plot(tmpsdata[zt],tmp[rev(zt)],col=cl[i],ylab=paste("moment = ",round(x@moment[i],2)), xlab="data",main=datanames[i],cex.lab=0.9)
			}
			mtext(paste("pmoments  : bootstrap"), side = 4, adj = 0, padj=0, col = "gray", cex = 0.4)
		}
		title(paste("Cumulative Sample Partial Moment-to-Threshold",sep=""), outer=TRUE, line=-1, cex=0.75)
	}
}

.manalyticalplot<-function(x,...)
{
	# sample versus analytical
	# 2 plots (density and sample with partial moment shading on threshold
	# and barplot with overlayed moments (sample vs analytical)
	
	print("multiple dataset/analytical method: no plot method for this type of data")
}