You can run this code directly at RStudio Cloud after a free registration and importing it into your own space. On the right panel named FILES, click on Fuego.Rmd to open it. Now you can tweak the code and run it by yourself but be patient the first time since a number of components will need to be installed. The perron.csv data source is included.

As a general rule, you can execute so-called “chunks” of R code by clicking the Run button within the chunk. The first time you open this file, you’ll need to install the following packages (remove the comment sign # before executing).

suppressMessages(library(tseries))
suppressMessages(library(forecast))
suppressMessages(library(FitAR))
suppressMessages(library(strucchange)) 
vcov.ri <- function(x, ...) kernHAC(x, 
kernel = "Quadratic Spectral", prewhite = 1, approx = "AR(1)", ...)
#setwd("~/Dropbox/Doks/Risque/Vegeta/data/")

First, we load the dataset, declaring its time series nature.

raw <- read.csv("perron.csv")
summary(raw)
##       year           risk        
##  Min.   :1961   Min.   : 0.2566  
##  1st Qu.:1976   1st Qu.: 1.4818  
##  Median :1991   Median : 2.8133  
##  Mean   :1991   Mean   : 3.8482  
##  3rd Qu.:2006   3rd Qu.: 4.1489  
##  Max.   :2021   Max.   :17.8117
DAT <- ts(raw, start=c(1961), end=c(2021),frequency=1)

The ARMA model selection based on the AIC criterion suggest 5 lags. We test for a unit root with Phillips-Perron. As we obtain a very low p-value, we may thus reject the null hypothesis that there is one unit root, meaning that forest fire risk is trend stationary (TS) over the entire period.

X = DAT[,"risk"]
plot(X)

Acf(X)

SelectModel(X, ARModel="AR", Criterion="AIC", Candidates = 15, Best=5)
##   p AIC-Exact AIC-Approx
## 1 5  152.0988  -7.825849
## 2 6  154.0679  -6.881227
## 3 4  154.5584  -9.792977
## 4 3  155.0112  -6.948193
## 5 7  155.1187  -5.588025
pp.test(X, alternative = c("stationary"))
## Warning in pp.test(X, alternative = c("stationary")): p-value smaller than
## printed p-value
## 
##  Phillips-Perron Unit Root Test
## 
## data:  X
## Dickey-Fuller Z(alpha) = -43.385, Truncation lag parameter = 3, p-value
## = 0.01
## alternative hypothesis: stationary

Since the forest risk time series displays a 5 year lag, we form the 5-year moving average that displays a clear break.

Y1 = rollmean(DAT,5,align="center")
Y = Y1[,"risk"]
plot(Y)

We test for a break in the Y trend using the first difference and applying a test for structural change. The Bayesian Information Criteria recommended by Bai & Perron suggests one break in 1980 with pre and post slopes of 0.46 and -0.22 (in per-thousand) though with limited significance (p-values of 2% and 10%).

U = diff(Y)
Z <- breakpoints(U ~ 1)
summary(Z) 
## 
##   Optimal (m+1)-segment partition: 
## 
## Call:
## breakpoints.formula(formula = U ~ 1)
## 
## Breakpoints at observation number:
##                      
## m = 1     17         
## m = 2   8 17         
## m = 3   8 17    34   
## m = 4   8 17 29 37   
## m = 5   8 16 24 34 42
## 
## Corresponding to breakdates:
##                                 
## m = 1        1980               
## m = 2   1971 1980               
## m = 3   1971 1980      1997     
## m = 4   1971 1980 1992 2000     
## m = 5   1971 1979 1987 1997 2005
## 
## Fit:
##                                              
## m   0      1      2      3      4      5     
## RSS  39.52  34.22  31.47  29.58  29.10  28.94
## BIC 147.46 147.44 150.79 155.39 162.52 170.27
plot(Z)

confid <- confint(Z,breaks = 1)
plot(U)
lines(Z)
lines(confid)
## Warning: Confidence intervals outside data time interval
##   from 1964 to 2019 (56 observations)

coef(Z, breaks = 1)
##             (Intercept)
## 1964 - 1980   0.4578288
## 1981 - 2019  -0.2114201
sapply(vcov(Z, breaks = 1, vcov = vcov.ri), sqrt)
## 1964 - 1980 1981 - 2019 
##   0.1776169   0.1098412
confint(Z, breaks = 1, vcov = vcov.ri)
## 
##   Confidence intervals for breakpoints
##   of optimal 2-segment partition: 
## 
## Call:
## confint.breakpointsfull(object = Z, breaks = 1, vcov. = vcov.ri)
## 
## Breakpoints at observation number:
##   2.5 % breakpoints 97.5 %
## 1     5          17     31
## 
## Corresponding to breakdates:
##   2.5 % breakpoints 97.5 %
## 1  1968        1980   1994
fac.ri <- breakfactor(Z, breaks = 1)
fm.ri <- lm(U ~ 0 + fac.ri)
plot(U)
lines(fitted(Z, breaks = 1), col = 2)

summary(fm.ri)
## 
## Call:
## lm(formula = U ~ 0 + fac.ri)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.04930 -0.32118  0.07086  0.30166  2.01096 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)  
## fac.risegment1   0.4578     0.1931   2.371   0.0213 *
## fac.risegment2  -0.2114     0.1275  -1.659   0.1030  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7961 on 54 degrees of freedom
## Multiple R-squared:  0.1343, Adjusted R-squared:  0.1022 
## F-statistic: 4.187 on 2 and 54 DF,  p-value: 0.0204