Introduction

This vignette illustrates a complete workflow for fitting theexponential factor copula model using theeFCM package,including data preprocessing, model fitting, diagnostic checking, andsimulation. We illustrate the workflow using precipitation data from acounterfactual climate scenario, chosen here purely as an exampledataset.

Exponential factor copula model

The exponential factor copula model (eFCM; Castro-Camilo and Huser,2020) is stochastically defined as

\[W(s) = Z(s) + V,\]

where:

• \(Z(s)\),\(s \in \boldsymbol{S} \subset\mathbb{R}^2\), is a standard Gaussian process with stationarycorrelation function\(\rho(h)\),
  
• \(h = \| s_1 - s_2 \|\) is thedistance between two locations\(s_1, s_2 \in\boldsymbol{S}\),
• \(V \sim \text{Exp}(\lambda)\) isan exponentially distributed random variable with rate parameter\(\lambda\), independent of both\(Z(s)\) and the spatial location\(s\).

The process W(s) belongs to the class of asymptotically dependentmodels and may be viewed as a Gaussian location mixture. Dependence isintroduced through the common factor V, which induces asymptoticdependence while allowing flexible representation of sub-asymptoticextremal dependence. Because this factor does not vary spatially, themodel is particularly well-suited to spatially homogeneous regions,where it is reasonable to assume similar marginal distributions and tailbehavior across sites.

If\(Z_j = Z(s_j)\), for\(j = 1,\ldots,d\), the multivariate latentGaussian vector\(\boldsymbol{Z} = (Z_1,\dots, Z_d)^\top\) follows a multivariate normal distribution,\(\boldsymbol{Z} \sim \Phi(\cdot;\Sigma_Z),\), where the covariance matrix\(\Sigma_Z\) is determined by the correlationfunction\(\rho(h)\). In thisapplication, we adopt the exponential correlation function,a specialcase of the Matérn class,\[\rho(h) = \exp\left(-\frac{h}{\delta}\right),\]

where\(\delta > 0\) is a rangeparameter controlling the spatial decay. With this, the jointdistribution of the process\(W\) is\[F_d^W(\mathbf{w}) =\Phi_{D}(\mathbf{w};\mathbf{\Sigma}_{\mathbf{Z}}) - \sum_{j =1}^D\exp\left(\frac{\lambda^2}{2} - \lambdaw_j\right)\Phi_{D}\left(\mathbf{q}_{j,0}-\mathbf{\mu}_{j,0};\mathbf\Omega_{j,0}\right),\] where\[\mathbf{q}_{j,0} = \begin{pmatrix}\mathbf{w}_{-j} -\mathbf{\Sigma}_{\mathbf{Z};-j,j}w_j\\ 0\end{pmatrix}, \quad\mathbf{\mu}_{j,0} = \begin{pmatrix} (w_j - \lambda)(\mathbf{1}_{d-1} -\mathbf{\Sigma}_{\mathbf{Z};-j,j})\\ \lambda -w_j\end{pmatrix},\]\[\mathbf\Omega_{j,0}=\begin{pmatrix} \mathbf{1}_{d-1}\mathbf{1}_{d-1}^T -2\mathbf{\Sigma}_{\mathbf{Z};-j,j}+\mathbf{\Sigma}_{\mathbf{Z};-j,-j}&\mathbf{\Sigma}_{\mathbf{Z};-j,j}-\mathbf{1}_{d-1}\\\mathbf{\Sigma}_{\mathbf{Z};-j,j}-\mathbf{1}_{d-1}&1\end{pmatrix},\] In particular, by setting\(d=1\), the marginal distribution may bewritten as\[F_1^{\mathbf{W}}(w;\lambda) =\Phi(w) - \exp(\lambda^2/2 - \lambda w)\Phi(w - \lambda).\]

Data Preparation

We start by loading the package and the pre-processed weeklycounterfactual precipitation data. Specifically,counterfactual contains weekly maxima of precipitationunder natural-only forcing, andLonLat provides spatialcoordinates for each station.

library(eFCM)# Load weekly precipitation maxima (pre-aggregated)data("counterfactual")# matrix/data frame: [weeks × stations]# Load station coordinatesdata("LonLat")coord<- LonLat

C:\2592713L3906b7fc4.R

We can visualize how spatially averaged weekly maxima evolve overtime:

plot(1:nrow(counterfactual),apply(counterfactual,1, mean),type ="l",xlab ="Week",ylab ="Mean Precipitation (mm)",main ="Weekly Maxima of Counterfactual Precipitation")abline(h =quantile(apply(counterfactual,1, mean),0.9),col ="red",lty =2)

C:\2592713L3906b7fc4.R

Create Spatial Data Objects

The functionfdata() converts the spatiotemporalprecipitation data into the format required for model fitting:

cf_data=fdata(counterfactual, coord,cellsize =c(1,1))

C:\2592713L3906b7fc4.R To reduce vignette build time, we loadprecomputed results:

data(cf_data)data(fit)

C:\2592713L3906b7fc4.R

Model Fitting

We fit the exponential factor copula model to the counterfactual datausingfcm():

fit=fcm(s =1, cf_data,boot = T,R =1000)

C:\2592713L3906b7fc4.R Here,s = 1 is the index of thegrid point and the and the exceedance threshold is set viathres, a quantile level in (0,1) (default 0.9), andboot = TRUE,R = 1000 requests 1,000 bootstrapreplications for uncertainty quantification.

Model Summary and Diagnostics

summary(fit)#>#> === Summary of Factor Copula Model Fit ===#> Grid location: ID = 1 | Longitude = -7.844 | Latitude = 37.120#>#> Number of neighbors: 2#> Neighbor coordinates:#>   - ID = 1 | Longitude = -8.438 | Latitude = 37.120#>#> Model Coefficients (95% CI):#> $`Estimate (Bootstrap)`#>          par      lower      upper#> 1   2.519319  0.5087895   35.00694#> 2 791.454198 13.5021393 1647.42038#>#>#> Objective value (negative log-likelihood):#> [1] -18.10454#>#> Information criteria:#>       AIC       BIC      AICc#> -32.20908 -24.29223 -32.17783#>#> Fitting arguments:#> $thres#> [1] 0.9#>#> $nu#> [1] 0.5#>#> $censorL#> [1] TRUE#>#> === End of Summary ===

C:\2592713L3906b7fc4.R

We can extract point estimates of the parameters. Recall that theparameter estimates (\(\lambda\) and\(\delta\)) describe the strength ofthe common factor and the spatial range of dependence, respectively:

coef(fit,method ="hessian")#>          par       lower      upper#> 1   2.519319   0.2104598   30.15764#> 2 791.454198 133.2100399 4702.34637coef(fit,method ="boot")#>          par      lower      upper#> 1   2.519319  0.5087895   35.00694#> 2 791.454198 13.5021393 1647.42038

C:\2592713L3906b7fc4.R

We can also compute the usual model selection criteria indices:

logLik(fit)#> [1] 18.10454AIC(fit)#>       AIC#> -32.20908BIC(fit)#>       BIC#> -24.29223AICc(fit)#>      AICc#> -32.17783

C:\2592713L3906b7fc4.R

Goodness-of-Fit

We can draw Q-Q plots to compare empirical and model-based upper tailquantiles for a given station. By default,qqplot() alsooverlays a generalized Pareto distribution (GPD) fit; this can besuppressed by settinggpd = FALSE.

qqplot(fit,which =1)

C:\2592713L3906b7fc4.R

We can also assess extremal dependence using the conditionalexceedance probability\(\chi_h(u)\),which measures the probability of simultaneous exceedances at high butfinite thresholds. For two locations\(s_1\) and\(s_2\) separated by distance\(h\), with respective vector components\(W(s_1)\) and\(W(s_2)\),\(\chi_h(u)\) is defined as\(\chi_h(u)=\lim_{u\to1}\Pr(W_{s_1}(W(s_1))>u\mid F_{s_2}(W(s_2))>u)\).The functionchiplot() ineFCM plotsmodel-based estimated of\(\chi_h(u)\)alongside their empirical counterpart for a range of values of\(u\). The package also provides two methodsfor pointwise confidence intervals: a normal approximation to the MLE orbootstrap resampling. Both approaches are illustrated below.

chiplot(fit,method ="hessian",ylim =c(0,1))

C:\2592713L3906b7fc4.R

chiplot(fit,method ="boot",ylim =c(0,1))

C:\2592713L3906b7fc4.R

Simulation from the Fitted Model

To generate synthetic precipitation fields consistent with theestimated extremal dependence structure, we can simulate from the fittedmodel as follows:

lbda<- fit$par[1]delta<- fit$par[2]dist<-rdist.earth(fit$coord)sim<-rmfcm(lbda, delta, dist,n =1e4)

C:\2592713L3906b7fc4.R

Summary

This vignette has demonstrated the complete pipeline for fitting theexponential factor copula model using theeFCM package. Itcovered preprocessing of data (in this case, climate model outputs),model estimation, diagnostic checks, and simulation. For a more detaileddiscussion of the method and its application to extreme eventattribution, see Li and Castro-Camilo (2025+).

References

Castro-Camilo, D., & Huser, R. (2020). Local likelihoodestimation of complex tail dependence structures, with application toU.S. precipitation extremes.Journal of the American StatisticalAssociation,115(531), 1037–1054.https://doi.org/10.1080/01621459.2019.1611584

Li, M., & Castro-Camilo, D. (2025+). On the importance of tailassumptions in climate extreme event attribution. arXiv preprint.https://doi.org/10.48550/arXiv.2507.14019