Dear All,
I have a question about package "ismev", its function "gpd.fit", and
interpretation of the results.
I used the package ismev to do an extreme value analysis on a fire
dataset. Two variables are used in the analysis. The focal variable is
acreage burned per fire, ranging from 1 to 5000 acres per fire. In
total, there are 69,980 observations. The date covers 1991 to 2008. I
created a time variable indexed by day, ranging from 1 (i.e.,
July/1/1991) to 6,708 (Nov/10/2008). In total, there 6,708 days. I used
function "gpd.fit" to estimate two threshold models, one without time
and the other with time (i.e., stationary vs. nonstationary). The key
outputs are attached as follows.
My first question is how to interpret the coefficient for time trend (0.
00473462). Does it mean the fire size increase by 0.00473462 acre per
day? My second question is related to calculating return levels. There
are multiple fires on individual calendar days. Each year there are 2000
to 5000 fires. I am not sure how to specify the number of observations
per year, i.e., the value of "npy" and which formulas should be used.
Thank you for your help.
Edwin
Changyou Sun, Ph.D.
Associate Professor
Natural Resource Policy & Economics
Box 9681, Department of Forestry
Mississippi State University
Mississippi State, MS 39762
#363 Thompson Hall
(662) 325 7271 (ph), 325 8726 (fax)
c...@cfr.msstate.edu
www.cfr.msstate.edu/forestry
> xa <- data$area; head(xa); NROW(xa)
[1] 1 2 1 1 2 1
[1] 69980
> ta <- as.matrix(data$time); head(ta); tail(ta); dim(ta)
[,1]
[1,] 1
[2,] 1
[3,] 1
[4,] 1
[5,] 1
[6,] 1
[,1]
[69975,] 6707
[69976,] 6707
[69977,] 6707
[69978,] 6707
[69979,] 6708
[69980,] 6708
[1] 69980 1
> th <- 100
> ra <- gpd.fit(xa, threshold=th) # stationary
$threshold
[1] 100
$nexc
[1] 1036
$conv
[1] 0
$nllh
[1] 5971.029
$mle
[1] 77.8220222 0.4094496
$rate
[1] 0.01480423
$se
[1] 3.71784848 0.03845162
> rb <- gpd.fit(xa, threshold=th, ydat=ta, sigl=1) # nonstationary
$threshold
[1] 100
$nexc
[1] 1036
$conv
[1] 0
$nllh
[1] 5966.575
$mle
[1] 61.60260272 0.00473462 0.40523571
$rate
[1] 0.01480423
$se
[1] 5.902059762 0.001556933 0.037975847
> lla <- -1*as.numeric(ra$nllh)
> llb <- -1*as.numeric(rb$nllh)
> dif <- 2*(llb-lla)
> lla; llb; dif
[1] -5971.029
[1] -5966.575
[1] 8.907359
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
