R help - Jan 2010 - Package ismev, gpd.fit, and interpretation for statistics of extreme values

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)

csun@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]]