Impact of Global Warming on Extreme Rainfall in Taiwan Cheng-Ching Lin Master Thesis National Tsing Hua University 2023 Abstract
The relationship between global warming and extreme rainfalls in Taiwan was examined in this study. Taiwan rainfall data from TCCIP, a project led by MOST, were analyzed. North
Hemisphere reference temperature data from NCEI led by NOAA. The yearly maximum of daily rainfall was focused on and the PGEV model, as proposed by Olafsdottir et al., was used to fit the extreme values and make inferences. The PGEV model integrates the General Extreme Value (GEV) and Peak over Threshold (PoT) approaches,
which are commonly used to analyze extreme data. Relative intensity and return value were used to show the connection between temperature and extreme rainfall. Results indicated that the intensity of extreme rainfall in Taiwan increases as the temper
ature rises. However, the effects of global warming on the frequency and intensity of extreme rainfalls varied by region. In the north and south regions, the frequency of extreme rainfalls
changed, while in the center and east regions, the intensity of extreme rainfalls changed. Fur
thermore, according to the return value analysis, extreme rainfalls are likely to occur more frequently in the future. To account for differences between locations, Gaussian Process was used to smooth the results obtained using the PGEV model. In addition, simulations using the Gaussian copula and Gaussian Process were conducted to determine the quantile confidence intervals for each
PGEV model. The simulations showed that all tests comparing with models with and without
covariates are significant. Keywords: Extreme Rainfall, Extreme Value Analysis, Gaussian Process, Climate change,
Global Warming
1 1 Introduction
Extreme value analysis draw more attention recently, especially in the fields of environment and finance. In finance, how to predict and prevent the damage of extreme financial risk is the main interest of risk management. Extreme value analysis is mainly used in quantifying tail distribution, finding the measurement of risks such as Value-at-Risk (VaR) and Expected
Shortfall (ES), and describing the tails properties in the markets of stocks,
Sheraz et al.), futures, and cryptocurrencies. Nolde & Zhou provided a comprehensive introduction to the application of extreme value analysis in finance. In environmental science, extreme value analysis is applied to characterize the occurrence and intensity of natural disasters caused by extreme weather, such as floods and earthquakes. Extreme value analysis on spatial data is a popular topic in environmental research. Blanchet et al. modeled the extreme snow depth in Switzerland and found the directional effect of snowfall. Davison
et al. introduced the classical modeling of spatial extremes, and Huser
& Wadsworth provided a more advanced introduction to spatial extremes analysis. There is two popular approaches in extreme value analysis. One is called the block maximum
approach, and it focuses on maximum (or minimum) data in each time blocking, such as the daily maximum rainfall. The other is called the threshold exceedance approach, and it focuses on the data higher than a certain threshold. For modeling, the block maximum approach uses the generalized extreme value (GEV) model to characterize the tail distribution and probability;
the threshold exceedance approach uses the Peak over Threshold (PoT) model which explores the intensity and frequency of the extreme events of interest.
More details on extreme value models refer to the textbooks and review papers,
which give an overall classical introduction to the concepts of extremes analysis. 1 According to the AR6 report from IPCC
(Intergovernmental Panel on Climate Change), the global temperature is getting higher every year, and the global temperature is already risen by 1 . 1 [◦] C in 2021 relative to the period 1850-1900. Being aware of such climate change, Olafsdottir et al. aims to answer the question that the frequency and intensity of extreme rainfall changes as temperature rises in the Northeastern United States. To use the annual daily maximum precipitations with consistent data quality, and to maintain the analysis ability of the threshold
exceedance approach, they combined GEV and PoT approaches with Poisson distribution to build anew model called PGEV to solve the problem. Using different settings in PGEV models,
they analyze the impact on extreme rainfall under various scenarios of global warming. With the data evidence, they found that, as the temperature rises, the intensity of extreme rainfall is fairly stable but extreme rainfall events happen more often in the Northeastern United States. Similar to Olafsdottir et al., this thesis try to answer the following question:
Whether is frequency or intensity of extreme rainfall events changes as the temperature becomes higher in Taiwan? To answer the question, this study analyzes Taiwan rainfall data sourced from the TCCIP
(Taiwan Climate Change Estimation Information and Adaptation Knowledge Platform) of the
MOST (Ministry of Science and Technology). The mission of the TCCIP project is to provide the climate change scientific data service to achieve the following goals, including strengthening the high-resolution model simulation ability in Taiwan, supporting local impact research, inte
grating climate adaptation services such as climate scenarios, risk information, and adaptation
tools, and developing integration of adaptation knowledge and implementation framework to expand the application possibilities of science-based solutions. Among the research in TCCIP,
two studies are relevant to extreme rainfall analysis in Taiwan. Tung et al. transformed daily precipitation into probability indices with cumulative distribution function
2 of GEV distribution and use probability indices to detect whether the probability of extreme rainfall events becomes higher when years go by. Henny et al. used TC
CIP gridded spatial rainfall data and defined five rainfall seasons to proceed the extreme value analysis under different rainfall types separately. They also defined “extreme rainfall events”
as the rainfall which is above the 99th percentile in each season and each location, and use such extreme rainfall data to fit the Theil-Sen slope. Their result reveals that, in Winter, Spring,
and Typhoon seasons, the frequency of extreme rainfall events in Taiwan becomes higher; the intensity of extreme rainfall events also becomes stronger in Typhoon season. This thesis adopts the PGEV approach to infer the relationship between temperature and the intensity and frequency of extreme rainfall in Taiwan. In par
ticular, the temperature is used as a temporal covariate. Based on the fitted PGEV model,
the analysis of relative frequency and return level of the extreme rainfall in Taiwan shall be investigated under some scenarios of temperature raise. In the end, the Gaussian process is used to smooth the PGEV parameters over the spatial region of interest for better visualization of parameter estimates and to quantify the uncertainty for modeling inference. This thesis is organized as follows: Section 2 introduces the temperature and rainfall data and the exploratory data analysis. The extreme value models are introduced in Section 3. The parameter estimation and detail fitting procedures are given in Section 4. Section 5 presents the application of extreme rainfall in Taiwan. Section 6 concludes the thesis contribution. The appendix provides the simulation and the confidence interval construction for some inferences in the application using the Gaussian process. 3 2 Data
This work uses temperature in NCEI and annual maxima of daily precipitation data in TCCIP
to process the following analysis. 2.1 Rainfall Data
The precipitation data in TCCIP is jointly created by the Meteorological Bureau, the Water
Resources Administration, the Civil Aviation Administration, the Environmental Protection
Administration, and the Forestry Research Institute, among other units. The daily rainfall data are available from 1960 to 2020, on a 5-kilometer resolution grid. There are 1311 5km ×
5km pixels covered Taiwan. The annual maximum of daily rainfall is extracted at each grid pixel and coupled with the latitude and longitude corresponding to the pixel. The annual maximum of daily precipitation in Taiwan during the period of 1960
2020 are shown in Fig 1. The TCCIP rainfall data are generated separately in 4 regions: North,
Center, South, and East. Since the boundary of these regions is overlapped, there are replicated measurements that occurred at boundary pixels (violet pixels) shown in left map of Fig 2.
To handle this problem, this study only keeps the maximum among replicated measurements at each boundary location, since the underlying analysis focuses on extreme rainfall. To explore data patterns, z t ( s j ) at each location s j is first viewed as a function of t . To visualize the temporal patterns of extreme rainfall events, the functional boxplot of the annual maximum daily rainfall among all pixels and that among pixels in each subregion is shown in Fig 3. The x-axis indicates years, and the y-axis indicates the rainfall measurements (mm). Here, the range of y-axes in all subplolts are. The black line in the middle represents the median function. The pink highlights the region of 25% to
75% quantiles among the functions at all pixels. The upper and lower blue curves indicate 1.5
times of pink area height. If at least one point in a curve is outside the upper or lower blue
4 1960 1970 1980 1990
2000 2010 2020
Rainfall (mm) 925
745
565
385
205
25
Figure 1: Annual maximum of daily rainfall in 5km × 5km resolution grid from TCCIP database in the period of 1960-2020. curve, the curve is detected as a potential outlier. From Fig 3, the median of annual maximum daily rainfall does not have clear decreasing or increasing trend patterns for the entire or local
regions. Generally speaking, the variation of annual maximum daily rainfall becomes larger and unstable in recent years. To see which region of Taiwan variates the most, Fig 3 shows the functional boxplot of different regions. In the South region, a high peak occurred in 2009,
due to typhoon Morakot, the deadliest typhoon in recorded history in Taiwan, hit the South region that year. In the Center region, the upper blue bound is far away from the center range
(pink region) in the period of year 2000-2010, which means most of the locations in the Center
5 (a) The map of regions and boundary pixels (b) The map shows the potential outliers pixels
Figure 2: The left map shows the regions and boundary pixels. The green, red, cyan, and orange pixels represent North, Center, South, and East regions respectively. The violet pixels represent the overlapping boundaries of 4 regions in Taiwan; The right map shows the potential outliers pixels detected in Fig 3. regions have a large variation than before. Finally, in the North and East regions, the trend and variation are quite stable than in the Center and South regions. There are three outlier candidates shown in red dashed lines. These outlier candidates are shown in right map of Fig
2 as red pixels. Two of candidates in Center region are at the boundary pixels, and the other one in East region is in the middle of Yilan. Finally, the temporal dependence in the series of annual maximum daily rainfall is checked. The sample ACFs are calculated for data at each pixel, and plotted as boxplot among all pixels for each lag in Fig 4. Since no box is situated outside the confidence band (shown
as the dotted lines), the annual maximum daily rainfall series can be assumed as temporal
uncorrelated. Namely, later in the modeling and analysis, the correlations between rainfall
6 Taiwan
1960 1970 1980 1990 2000 2010 2020
North
1960 1970 1980 1990 2000 2010 2020
Center
1960 1970 1980 1990 2000 2010 2020
South
1960 1970 1980 1990 2000 2010 2020
East
1960 1970 1980 1990 2000 2010 2020 Year
Figure 3: The functional boxplot of different regions of Taiwan’s and whole Taiwan regions annual daily maximum rainfall in the period of 1960-2020. 7 1 2 3 4 5 6 7 8 9 11 13 15 17 19
Lag
Figure 4: The x-axis means time lag, the y-axis means ACF of rainfall, and a box means the distribution of ACF at a certain lag in all pixels in Taiwan. For each lag, there are 1311 ACF values to construct a box.
observations among different years can be ignored. 2.2 Temperature Data
The thesis uses the North Hemisphere annual average reference temperature data provided by the National Centers for Environmental Information (NCEI), supported by the National
Oceanic and Atmospheric Administration (NOAA). NCEI is the leading authority on environ
mental data. They maintain one of the largest archives of environmental research around the
world. The temperature data refers to the average temperature during 1901-2000. To explore the long-term trend patterns, the yearly reference temperature are further smoothed based on the Lowess method), shown as the red line in Figure 5. This “smoothed temperature” is considered as a
covariate to study the effect of climate change in extreme rainfall analysis. 8 1880 1900 1920 1940 1960 1980 2000 2020
Year
Figure 5: For each point is the annual average reference temperature in north hemisphere with respect to 1901-2000. The x-axis is years, and the y-axis is the reference temperature( [◦] C).
The red line is the smoothed reference temperature trend using Lowess method. 2.3 Notations
Finally, define the notations for data:
- each pixel in the grid is represented as j
- latitude and longitude of each pixel are represented as a coordinate vector s j
- t represents year
- annual maximum daily rainfall at location s j and year t is written as Z t ( s j )
- the vector of observed yearly maximum of daily rainfall at location s j from year 1 to T
is z ( s j ) ≡ , . . ., z T ( s j )) [′]
- the smoothed reference temperature at year t is denoted as x t
9 3 Models for Extreme Value Analysis
In this section, two fundamental models in extreme value analysis, GEV and PoT, are intro
duced. The PGEV model, which integrates the formulation of both
GEV and PoT models, will be introduced afterward. 3.1 Generalized Extreme Value Models (GEV)
For the GEV model, the data it fits are block maximum data. The maximum statistic among the observations in a block of time, such as month or year, follows the GEV model asymptotically. The CDF of GEV( µ, σ, γ ) distribution satisfies
G ( z ; µ, σ, γ ) = exp
 −
�
1 + γ [z][ −] [µ]
σ
− 1 /γ []
, γ ̸ = 0 ,
�

for z satisfying γ ( z − µ ) < σ, where µ ∈ R is the location parameter, σ > 0 is the scale
parameter, and γ ∈ R is the shape parameter. The corresponding density of GEV( µ, σ, γ ) is f G ( z ; µ, σ, γ ) =
σ
for z satisfying γ ( z − µ ) < σ . �
− 1 /γ []
,
�

1 + γ [z][ −] [µ]
 −
�
1 + γ [z][ −] [µ]
σ
− 1 /γ− 1
exp �
σ
If the shape parameters γ = 0, the GEV model becomes a Gumbel model:
��
, z ∈ R . G ( z ; µ, σ ) = exp
�
− exp
�
− [z][ −] [µ]
σ
The parameters in the GEV model have a close connection to those in the PoT model. Their connections are given in Section 3.3. 3.2 Peak Over Threshold Models (PoT)
For the PoT model, the data it fit are the exceedance
values over a high threshold c . The occurance of such over-threshold events follows a Poisson
10 process with the rate parameter λ c depending on c . The exceedances are assumed to be mu
tually independent, independent of the Poisson process, and follow the PoT( σ c, γ ) model. In
particular, the exceedance ( Z d − c ) follows the generalized Pareto distribution with the CDF:
H d ( z ; σ c , γ ) ≡ P ( Z d < c + z|Z d > c ) = 1 −
�
1 + γ [z]
σ c
− 1 /γ
, σ c + γz > 0 , �
where σ c > 0 is a scale parameter depending on the threshold c, and γ is a shape parameter. The corresponding density of H d ( z ; σ c , γ ) is
− 1 /γ− 1
, σ c + γz > 0 . �
f H ( z ; σ c , γ ) =
σ c
�
σ c
1 + γ [z]
If the shape parameter γ = 0, the PoT model becomes an exponential model:
�
�
, z > 0 . H d ( z ; σ c ) = 1 − exp
− [z]
σ c
3.3 PGEV Models
This work aims to analyze the frequency and intensity of extreme rainfall events with annual daily maximum rainfall, but it cannot be done by the GEV model alone. Since the frequency and intensity of extreme rainfall events can be represented respectively as the rate of rainfall events above thresholds and the exceedences values from thresholds, these inference properties are came from PoT model not GEV model. Therefore, Olafsdottir et al. introduced the PGEV model, which fits the block maximum data but infers the parameters in
PoT models. The logic is illustrated as follows. Suppose the extreme daily rainfall over a threshold c followed by a Poisson process with a yearly rate λ c .
Let N be the random number of excesses in a year and let Z be the annual
11 maximum of daily rainfall. It follows that
P ( Z < z ) =
≈
=
365 � P ( Z < z|N = n ) P ( N = n )
n =0
∞ � [ H d ( z − c ; σ c , γ )] [n] [ λ] n c [n] ! [e] [−][λ] [c]
n =0
∞ �
n =0
1 + γ [z][ −] [c]
n

1 n ! 
 [λ] [c]
1 −
�

 

e [−][λ] [c]
σ c
− 1 /γ [] �

1 −

 

 [e] [−][λ] [c]
= exp
= exp
= exp
= exp
= exp
 
 [λ] [c]
      
     

1 + γ 
�
1 + γ [z][ −] [c]
σ c λ ~~[γ]~~ c
σ c
− 1 /γ [] �

 
 [−][λ] [c]
�
1 + γ [z][ −] [c]
σ c
− /γ []  �

 
 [−]
 
 [−]
�
�
1 + λ [−] c [γ] + γ [z][ −] [c]
− 1 /γ []  �
− 1 σ c

γ 1 + λ ~~[γ]~~ c [σ] c
�
[c] c [σ] [c]
γ [−] [λ] [γ] γ
λ [γ] c [(] [z][ −] [c] [) +] [ σ] [c]
γ
− 1 /γ []  ��

− 1 /γ []  ��
− 1 /γ []
     
,

     

 
z −
�
λ [γ] c [−] c + σ c
γ
�
−
in which the third equation is due to. Olafsdottir et al. called the CDF of a PGEV( σ c, γ c ) model given the threshold c :
− 1 /γ []
     
.      

z −
12
�

 ~~~~ 
λ [γ] c [−] c + σ c
�
F c ( z ) = exp
      
     

1 + γ 
σ c λ ~~[γ]~~ c
γ
− The corresponding density of PGEV is
    
 − 1 /γ− 1
 ~~~~ 
exp


 ~~~~ 
− 1 /γ []     
1 f c ( z ) = σ c λ ~~[γ]~~ c

λ [γ] c [−]
 z − c + σ c
� γ 1 + γ
σ c λ ~~[γ]~~ c
 
σ c λ ~~[γ]~~ c
γ
�
�
    
−

    
λ [γ] c [−]
 z − c + σ c
� γ 1 + γ
σ c λ ~~[γ]~~ c
 
σ c λ ~~[γ]~~ c
γ
. Compare with, PGEV is a special case of GEV, i.e., GEV( µ GEV, σ GEV, γ ), but
parametrized via ( λ c , σ c , γ ) given the threshold c, satisfying
λ [γ] c [−] µ GEV = c + σ c ~~,~~ σ GEV = σ c λ [γ] c [.]
γ
To be more general, this thesis allows the covariate effect in the rate and scale parameters and considers the regression form for ( λ c , σ c ) at each pixel j and year t as:
log λ c ( s j , x t ) = β 0 j + β 1 j x t ,
log σ c ( s j , x t ) = α 0 j + α 1 j x t ,
where x t is the long-term temperature given in Section 2.2 shown in Fig 5. Consequently, are allowed to be pixel-wise and year-specific forms:
µ GEV ( s j , x t ) = c j + σ c ( s j , x t ) [λ] [c] [(] [s] [j] [, x] [t] [)] [γ] [j] [ −] ~~,~~ σ GEV ( s j , x t ) = σ c ( s j , x t ) λ c ( s j , x t ) [γ] [j] . γ j
3.4 Temperature Effects on Intensity of Rainfall
This work evaluates the relative frequency change, relative scale change and quantifies the return level subject to temperature raises to investigate the impact of climite change on the distribution of the annual maximum daily rainfall in Taiwan. The ideas about the relative frequency change and the return level quentification are adopted from the approach by Olafsdottir et al.. 3.4.1 Relative Frequency and Scale Changes
According to and, the temperature effects the PGEV parameter ( λ c , σ c ). Let ∆ x =
x high − x low be the temperature change. The relative changes in the frequency parameter and
13 the scale parameter at the pixel j subject to temperature change ∆ x are
∆ λ ( s j ) ≡ [λ] [c] [(] [s] [j] [, x] [high] [)] [ −] [λ] [c] [(] [s] [j] [, x] [low] [)] = exp − 1 ,
λ c ( s j , x low )
∆ σ ( s j ) ≡ [σ] [c] [(] [s] [j] [, x] [high] [)] [ −] [σ] [c] [(] [s] [j] [, x] [low] [)] = exp − 1 .
σ c ( s j , x low )
Different scenarios of ∆ x will be discussed in the application. 3.4.2 Return Levels
A return level (denoted as R q ) of the annual maximum daily rainfall with a return period of
1 /q years is defined as a high threshold whose probability of exceedance is q ∈, i.e., the
-th quantile of the underlying extreme distribution. For example, for q = 0 . 05, R 0 . 05
means the intensity could occurred once in a 20-year period. For the GEV( µ, σ, γ ) model,
the return level satisfies
R q ≡ G [−] = µ − [σ]
γ
1 − [−][γ]
� �
, 0 < q < 1 . Plugging in the parameter estimates in Eq, the return level under the temperature x t at pixel j becomes
R q ( s j , x t ) = µ GEV ( s j , x t ) − [σ] [GEV] [(] [s] [j] [, x] [t] [)] { 1 − [−][γ] [j] }
γ j
�
− [λ] [c] [(] [s] [j] [, x] [t] [)] [γ] [j] [σ] [c] [(] [s] [j] [, x] [t] [)] { 1 − [−][γ] [j] }. γ j
=
�
c j + [σ] [c] [(] [s] [j] [, x] [t] [)]
~~(~~ λ c ( s j , x t ) [γ] [j] − 1) γ j
To evaluate how the temperature raise affects the intensity of extreme events, the proba
bility of a storm at pixel j having the rainfall intensity exceeded over R q ( s j , x low ) subject to a temperature change ∆ x = x high − x low is defined as
P [ Z t ( s j ) > R q ( s j , x low ) |x t = x high ]
− 1 /γ j []
     
λ c ( s j , x high ) [γ] [j] σ c ( s j , x high )
14
     
�
R q ( s j , x low ) −
�
c j + [σ] [c] [(] [s] [j] [, x] [high] [)]
[ λ c ( s j , x high ) [γ] [j] − 1] γ j

 ~~~~ 
= 1 − exp
     
−

     

1 + γ j If covariate has no effects on parameters, i.e., β 1 j = α 1 j = 0, Eq reduces to q . If the
covariate only effects on the rate parameter, i.e., α 1 j = 0 , β 1 j ̸ = 0, Eq is simplified as
�
= 1 − [∆] [λ] [(] [s] [j] . P [ Z t ( s j ) > R q ( s j , x low ) |x t = x high ] = 1 − exp
λ c ( s j , x high )
λ c ( s j , x low )
�
λ c ( s j , x low ) [ −] [q] [)]
In particular, Eq only relies on the relative changes in frequency parameter λ defined in Eq. On the other hand, if the covariate only effects on the scale parameter, i.e.,
α 1 j ̸ = 0 , β 1 j = 0, Eq is simplified as
P [ Z t ( s j ) > R q ( s j , x low ) |x t = x high ]
− 1 /γ j []
 �
 [.]
= 1 − exp
 
 [−] [exp]
�
1 −
�1 − exp) [−][γ] [j] [�]
∆ σ ( s j ) + 1
Eq means that if the relative intensity changes by ∆ σ ( s j ), for temperature ∆ x increase,
then a R q -level storm will happen with Eq. That is, the probability of R q -level storms in x low increase under x high . 15 4 Estimation
This study assumes the nonstationary PGEV model with spatial and temporal varying param
eters ( λ c ( s j , x t ) , σ c ( s j , x t ) , γ j ) at year t and pixel j satisfying and. For simplicity, this work only considers model fitting at each pixel using the maximum likelihood estimation. 4.1 Determining Threshold
To fit the PGEV model, the threshold c j for daily rainfall of each pixel j is needed to be determined in in advance. For simplicity, this thesis determines the thresholds c j by fitting
GEV with extreme value rainfall for each pixel j and ignoring the effect of temperature.
That
is, find the MLE of parameters from pdf of GEV Eq. The threshold is determined as follows
that, an extreme rainfall event is defined as rainfall that only happens on average in about
365 . 25 days. The PGEV parameters in Eq can be transformed into the formula
below:
−γ ˆ j p ) c j = ˆ µ j − [σ] [ˆ] [j] [ −] ˆ [λ] ~~,~~ γ j where λ p = 365 . 25. 4.2 Maximum Likelihood Estimation for PGEV
The study use maximum likelihood (ML) to estimate the parameters θ j = [′]
in PGEV model for each pixel j . Suppose the last temporal index is T . According to, the
16 log-likelihood of θ j given the observations z ( s j ) ≡ , .
. ., z T ( s j )) [′] at the j th pixel is
ℓ ( θ j ) =
= −
−
+
T �
t =1
T �
t =1
T �
t =1
T � log f c j ( z t ( s j ))
t =1
log
�
��
1 + γ j � �
1 +
γ j
�
exp [−] [/γ] [j] [�] �
. The ML estimate of θ j is defined as
ˆ θ j = arg max ℓ ( θ j ) ,
which is found using the maxLik package in R. 4.3 Hypothesis Testing and Model Selection
To find out whether the frequency or the intensity of extreme rainfall changes with temperature,
this work fits four PGEV models as Eq at each pixel j and finds MLE of parameters by solving log-likelihood Eq. The four models are:
- PGEV λ,σ model: the frequency and scale parameters varying with temperature and the parameter space is
Ω λ,σ = { θ : β 0 ∈ R , β 1 ∈ R , α 0 ∈ R , α 1 ∈ R , γ ∈ R } . The correspondence MLE is denoted as θ [ˆ] j [(] [λ,σ] [)] = arg max ℓ ( θ j ) θ j ∈ Ω
- PGEV λ model: only the frequency parameter changes with temperature and the param
eter space is
Ω λ = Ω ∩{α 1 = 0 } . The correspondence MLE is θ [ˆ] j [(] [λ] [)] = arg max ℓ ( θ j ) . θ j ∈ Ω λ
17 - PGEV σ model: only the scale parameter changes with temperature and the parameter space is
Ω σ = Ω ∩{β 1 = 0 } . The correspondence MLE is θ [ˆ] j [(] [σ] [)] = arg max ℓ ( θ j ). θ j ∈ Ω σ
- PGEV 0 model: the parameters is time-invariant and the parameter space is
Ω 0 = Ω ∩{β 1 = α 1 = 0 } .
The correspondence MLE is θ [ˆ] j = arg max ℓ ( θ j ). θ j ∈ Ω 0
To select the appropriate model, this study conducts three Likelihood Ratio tests (LRT)
and gets those p-values for each pixel j . The three tests are:
Test 1 H 0 : θ j ∈ Ω 0 vs H 1 : θ j ∈ Ω λ,
Test 2 H 0 : θ j ∈ Ω 0 vs H 1 : θ j ∈ Ω σ,
Test 3 H 0 : θ j ∈ Ω 0 vs H 1 : θ j ∈ Ω λ,σ . Accordingly, the LRT test statistics for Tests 1-3 and their null distribution under H 0 are given as follows respectively:
Test 1
Λ λ = − 2 ℓ] [)] [ −] [ℓ] [(ˆ] [θ] j [(] [λ] [)] [)] ∼ χ 1 [,] � �
Test 2
Λ σ = − 2 ℓ] [)] [ −] [ℓ] [(ˆ] [θ] j [(] [σ] [)] [)] ∼ χ 1 [,] � �
Test 3
Λ λ,σ = − 2 ℓ] [)] [ −] [ℓ] [(ˆ] [θ] j [(] [λ,σ] [)] ) ∼ χ 2 [.] � �
18 Besides LRTs, AIC) is also used for model
selection. Let θ [ˆ] be a generic notation for MLE of one of PGEV models. The definition of AIC
is:
AIC ( θ [ˆ] ) = − 2 ℓ ( θ [ˆ] ) + 2 dim( θ [ˆ] ) ,
where dim( θ [ˆ] ) denotes the number of parameters in θ [ˆ] . AIC will be calculated for the models fitted to each pixel. The lower value of AIC, the better fitting of models.
The whole fitting procedure is outlined in Algorithm 1. Algorithm 1 Steps for finding MLE of PGEV
Require: p for determining thresholds c j ;
1: for each j ∈ do 2: Fit GEV( µ, σ, γ ) with annual maximum data z ( s j ) to get ˆ µ j,ˆ σ j,ˆ γ j . 3: Substitute the parameters from previous step into Eq, and determine the threshold c j . 4: Set ˆ µ j,ˆ σ j,ˆ γ j from step 2 as initial value parameters of PGEV model, and use BFGS and Eq to find the MLE θ [ˆ] j . 5: For all models with covariates x t, perform LRT on the stationary models, get the p-value for each model and calculate AIC for each model. 6: end for
19 µ
361
334.2
307.4
280.6
253.8
227
200.2
173.4
146.6
119.8
93
(a) location parameters ˆ µ j
σ
161
148.3
135.6
122.9
110.2
97.5
84.8
72.1
59.4
46.7
34
(b) scale parameters ˆ σ j
γ
0.4
0.2
0.0
−0.2
−0.4
(c) shape parameters ˆ γ j
Figure 6: Parameter estimates fitted by GEV. 5 Application
5.1 GEV Parameters and Threshold
This work uses the annual maximum daily rainfall data z t ( s j ) from TCCIP to analysis the extreme rainfall property in Taiwan. The temperature data x t is from NCEI. The period of data is 61 years from 1960 to 2020 and j = 1 , 2 , ..., 1311 in a 5km x 5km grid. First, fit GEV( µ j, σ j, γ j ) models as Eq ignoring the effect of temperature for each pixel j . The estimated parameters of GEV are shown in Fig 6 as maps.
From Fig 6, the left map shows the location parameters µ j on each pixel in Taiwan, the middle of the South area has larger values, indicating those areas have a larger yearly maximum of daily rainfalls; the middle map shows the scale parameters σ j, the South and the East regions have larger values, showing the variation in these two regions are higher than other regions; the right map shows the shape parameters γ j, have a clear difference between Eastern and Western of Taiwan. After getting the parameters from GEV models, the thresholds c j for PGEV models are determined using Eq. Here, the thresholds are determined by the quantile of p = 0 . 99 in a GEV stationary model. The distribution of thresholds among pixels is shown in Fig 7. The histogram of thresholds is right-skewed, that represents few places having large extreme rainfall
20 50 100 150
Threshold
(a) The histogram of thresholds
162
147.8
133.6
119.4
105.2
91
76.8
62.6
48.4
34.2
20
(b) The map of thresholds
Figure 7: The histogram of thresholds is the left one, the dashed line is the median of thresholds. The right one is the map of thresholds. only. According to the map on the right of Fig 7, the southeast area and Yilan have large
thresholds, and this pattern is consistent with the rainfall maps from Fig 1. 5.2 PGEV Fitting and Model Comparison
After determining thresholds c j, the four PGEV models are fitted with the data z ( s j ) for each pixel j and the best model is determined by the test results of LRT Λ λ, Λ σ, and Λ defined in Eq
– Eq respectively. These tests are conducted for all pixels in Taiwan. The p-values of each test are collected. For each test, the p-values among pixels are sorted in ascending order,
denoted as p ( i ), and plotted against i/ 1311 as a curve in Fig 8. If one of the three tests is significant in most of the pixels in Taiwan, the representative curve will be concave-down and vice versa. From Fig 8, PGEV σ is most significant compared to the other models since the curve of Λ σ is most concave-down. The confidence band of each QQ-plot curve is constructed via the bootstrap method. The detailed procedure is provided in Appendix II. The curves with
21 1.00
0.75
0.50
Λ σ
Λ λ
Λ λσ
0.25
0.00
0.00 0.25 0.50 0.75 1.00
Ordered index
Figure 8: This is the QQplot of the p-value from LRT. The x-axis is ordered index and normalize it into. The black line is diagonal line which means not significant at all. The colored line means LRT p-value of each test. The gray line is 0.05. the corresponding confidence band are shown in Fig 9. From top left of Fig 9, all lines are inside the confidence band.
According to the confidence band, the p-values from three tests are not significantly superior to others. Moreover, the best model can be selected by AIC defined in Eq. Since the original data come from 4 regions: North, Center, South, and East, the fitting performance of 4 fitted PGEV
models in each region is summarized in Table 1. Also, Table 2 put the PGEV 0 models into the
comparison. The frequency given in the table is the counts among total pixels in each subregion
having the smallest AIC values among 4 fitted PGEV models. From Table 1, PGEV σ model is the most favored model with the lowest AIC in the Center and East regions, while PGEV λ
model is in the North and South regions. From Table 2, PGEV 0 is most favored in all regions except south. In South region, PGEV λ is most favored model. However, PGEV σ is favored in more pixels than PGEV λ . Fig 10 shows the pixels with lowest AIC models. Lots of pixels in
22 1.00
0.75
0.50
0.25
0.00
1.00
0.75
0.50
0.25
0.00
0.00 0.25 0.50 0.75 1.00
Ordered index
Λ λ
0.00 0.25 0.50 0.75 1.00
ordered index
Λ σ Λ λ Λ λσ
Λ σ
1.00
0.75
0.50
0.25
0.00
1.00
0.75
0.50
0.25
0.00
0.00 0.25 0.50 0.75 1.00
ordered index
Λ λσ
0.00 0.25 0.50 0.75 1.00
ordered index
Figure 9: The simulation result is the band in the picture above, and the previous lines are the same as seen in Fig 8. The top left is the QQplot of p-values, and the others are the QQplot of p-values with confidence band for each test. 23 Table 1: The frequency of the smallest AIC PGEV models in different regions
Regions PGEV λ PGEV σ PGEV λ,σ Total pixels North 141 99 15 255
Center 141 281 12 434
South 176 117 11 304
East 101 208 9 318
Total 559 705 47 1311
Table 2: The frequency of the smallest AIC PGEV models in different regions with PGEV 0
Regions PGEV 0 PGEV λ PGEV σ PGEV λ,σ Total pixels North 198 34 18 5 255
Center 249 47 128 10 434
South 15 171 107 11 304
East 200 26 85 7 318
Total 662 278 338 33 1311
the south regions are favored PGEV λ, while PGEV σ is favored in Yilan and the center region. Moreover, those pixels favored PGEV σ are clustered as band across Taiwan from west to east. Similarly, QQ plot on p-values for each subregion is shown in Fig 11.
The p-value in the south is particularly significant compared to those in the other regions, while the north is the least significant region. The results show that extreme rainfalls are more likely to be affected by
temperature, no matter freqency or intensity, in the South region. 5.3 Inference on Climate Change
Since Fig 11 and Table 1 shows that PGEV σ is significant in Center and East regions, while
PGEV λ is significant in North and South regions, this thesis infer the property of extreme rainfall in each model. 5.3.1 Inference Under PGEV σ Model
Since PGEV σ is picked as the best model in Center and East regions, the main focus of this thesis is on the parameter ˆ α 1 j, which represents the impact of temperature on the intensity of
24 PGEV 0 PGEV λ PGEV σ PGEV λσ
Figure 10: This is the map of pixels with lowest AIC models. extreme rainfall events. To see the distribution of ˆ α 1 j, Fig 12 shows the histogram of ˆ α 1 j for the entire Taiwan region and in 4 subregions. From Fig 12, the median of ˆ α 1 j (denoted as the solid red line) is positive in all regions. Also, more than half of the pixels are experiencing an increase in the intensity of extreme rainfall subject to temperature raise. Moreover, all pixels in the south have positive ˆ α 1 j, indicating a substantial impact of temperature raise on extreme rainfall in this region. In contrast, the median in the north is near to 0, showing less evidence of the climate impact on extreme rainfall events. To bring the location information into the PGEV σ model, and also to make plots spa
tially smoothed, a Gaussian Process (GP) is post-assumed to α 1 ( s j ) ≡ α ˆ 1 j to perform a high resolution map of α 1 ( s ) predict by 1km × 1km gridded pixels in Taiwan using kriging, displayed in Fig 13. The detailed process of constructing the GP model is provided in Appendix I. From
Fig 13, α 1 ( s ) is larger in the southern. It means that the temperature has more impact on extreme rainfall events once the temperature becomes higher. In contrast, the Yilan area has smaller value of α 1 ( s ), indicating less impact of temperature on extreme rainfalls events in this
25 Λ σ Λ λ Λ λσ
1.00
0.75
0.50
0.25
0.00
1.00
0.75
0.50
0.25
0.00
North
0.00 0.25 0.50 0.75 1.00
Ordered index
South
0.00 0.25 0.50 0.75 1.00
Ordered index
1.00
0.75
0.50
0.25
0.00
1.00
0.75
0.50
0.25
0.00
Center
0.00 0.25 0.50 0.75 1.00
Ordered index
East
0.00 0.25 0.50 0.75 1.00
Ordered index
Figure 11: This is the QQplot of the p-value from LRT but separate into 4 regions. 26 Taiwan
−1.5 −1.0 −0.5 0.0 0.5 1.0 1.5
North
−1.5 −1.0 −0.5 0.0 0.5 1.0 1.5
Center
−1.5 −1.0 −0.5 0.0 0.5 1.0 1.5
South
−1.5 −1.0 −0.5 0.0 0.5 1.0 1.5
East
−1.5 −1.0 −0.5 0.0 0.5 1.0 1.5
α [^] 1
Figure 12: The histogram of ˆ α 1 j in PGEV σ models in different regions. The dashed line means 0, and the solid line means the median of this distribution. 27 1.0 0.5 0.0 −0.5 −1.0
(a) The Taiwan map of ˆ α 1 j in PGEV σ models
2
1
0
−1
−2
(b) The Taiwan map of t-value of ˆ α 1 j in PGEV σ models
Figure 13: Left plot is spatially smoothed map of α 1 ( s ). Right plot is the t-value of α 1 ( s ) in each pixel.
According to AR6, the temperature changes are 1.5,
2, 3, and 4 [◦] C in near future relative to the period 1850-1900. Therefore, the temperature raises considered here are ∆ x = 0 . 5 , 1 , 2, and 3 [◦] C for the following analysis. Fig 14 shows the distribution of ∆ σ among all pixels as a function of ∆ x . The relative intensity ∆ σ is defined in Eq. The dotted line is the median of the individual distributions. Each line represents a temperature change scenario. When ∆ x rises, the median of relative intensity ∆ σ goes higher. For example, the median of ∆ σ is around 16% as ∆ x = 0 .
5 [◦] C, and the median of ∆ σ is around 140% as ∆ x = 3 [◦] C. Fig 15 shows the probability of extreme storm variation defined in Eq for q = 0 . 05 under
4 temperature change scenarios. From Fig 15, the distribution of 5% extreme storms becomes
28 Taiwan
1.5
1.0
0.5
Increase Temperature
0.5
1
2
3
0.0
0 2 4 Relative Frequency Increase
Figure 14: Each curve is a distribution of relative intensity increase for each pixels under certain situation that temperature changes. The dashed line is the median of each curve. The x-axis is relative frequency change, and the y-axis is density of the distribution. more often to happen as ∆ x increases. That is, if the temperature continues to increase in the future, the probability of happening extreme rainfall will become more than 5%. In other
words, originally the extreme rainfalls which happened every 20 years on average will occur more frequently when the temperature goes higher. It means in the future, the center and east of Taiwan will more likely to face extreme rainfalls than before, and the number of disasters followed by extreme rainfalls will also be increased, such as landslides and floods. 5.3.2 Inference Under PGEV λ Model
Since PGEV λ is picked as the best model in North and South regions, the main focus of this thesis is on the parameter β [ˆ] 1 j, which represents the impact of temperature on the frequency of extreme rainfall events. To see the distribution of β [ˆ] 1 j, Fig 16 shows the histogram of β [ˆ] 1 j for the entire Taiwan region and in 4 subregions. From Fig 16, the median of β [ˆ] 1 j (denoted as the solid red line) is positive in center and south regions. Moreover, most pixels in the south have positive β [ˆ] 1 j, indicating a substantial impact of temperature raise on extreme rainfall in
29 1.00
0.75
0.50
0.25
0.00
0.5 1 2 3 Relative Temperature Increase
Figure 15: The plot shows how much temperature rises, and the probability of extreme storm variation. The x-axis is the scenario of temperature increase, the y-axis is the probability of a 5% storm under a certain scenario of temperature increase, and each whisker box is a distribution of probability that extreme storm happens when temperature changes, the dashed line is 5% probability. this region. In contrast, the median in the north and east are near 0, showing less evidence of the climate impact on extreme rainfall events. The same as Fig 13 map, a Gaussian Process (GP) is used to perform a high-resolution map of β ( s ) in Fig 17. From Fig 17, β ( s ) is larger in the southern. It means that the temperature has more impact on extreme rainfall events once the temperature becomes higher. In contrast,
the Hualien area has smaller value of β ( s ), indicating less impact of temperature on extreme rainfalls events in this region. Fig 18 shows the distribution of ∆ λ among all pixels as a function of ∆ x . The relative frequency ∆ λ is defined in Eq. The dotted line is the median of the individual distributions. Each line represents a temperature change scenario. When ∆ x rises, the median of relative frequency ∆ λ goes higher. For example, the median of ∆ λ is around 12% as ∆ x = 0 .
The results is the same as Fig 15. 5.4 Hypothesis Testing on Models with Covariate
Since Fig 9 shows that the confidence bands of Λ λ, Λ σ, and Λ overlap each other, there are two more LRTs to confirm the relation between PGEV λ, PGEV σ, and PGEV λ,σ . The two LRT for each pixel j are:
Test a H 0 : θ j ∈ Ω λ vs H 1 : θ j ∈ Ω,
Test b H 0 : θ j ∈ Ω σ vs H 1 : θ j ∈ Ω. 32 Taiwan
1.5
1.0
0.5
Increase Temperature
0.5
1
2
3
0.0
0 2 4 Relative Frequency Increase
Figure 18: Each curve is a distribution of relative frequency increase for each pixels under certain situation that temperature changes. The dashed line is the median of each curve. The x-axis is relative frequency change, and the y-axis is density of the distribution. Accordingly, the LRT test statistics for Tests a-b and their null distribution under H 0 are given as follows respectively:
Test a
˜Λ λ = − 2 ℓ ( θ [ˆ] j [(] [λ] [)] [)] [ −] [ℓ] [(ˆ] [θ] [j] [)] ∼ χ 1 [,] � �
Test b
˜Λ σ = − 2 ℓ ( θ [ˆ] j [(] [σ] [)] [)] [ −] [ℓ] [(ˆ] [θ] [j] [)] ∼ χ 1 [.] � �
The results of LRTs from Test a-b and corresponding confidence bands are shown in Fig 20. The curve represents the p-value of each test in each pixel j . All lines are close to the diagonal black line, and it means PGEV λ and PGEV σ are not significant compared with PGEV λ,σ,
especially all confidence bands include the diagonal line. The results can be told from the confidence bands including the black diagonal line.
Therefore, we can conclude that the models with covariates are not significant to each other. 33 1.00
0.75
0.50
0.25
0.00
0.5 1 2 3 Relative Temperature Increase
Figure 19: The plot shows how much temperature rises, and the probability of extreme storm variation. The x-axis is the scenario of temperature increase, the y-axis is the probability of a 5% storm under a certain scenario of temperature increase, and each whisker box is a distribution of probability that extreme storm happens when temperature changes, the dashed line is 5% probability. 34 1.00
0.75
0.50
0.25
0.00
1.00
0.75
0.50
0.25
0.00
Λ [~] σ
0.00 0.25 0.50 0.75 1.00
ordered index
Λ [~] σ Λ [~] λ
Simulation
1.00
0.75
0.50
0.25
0.00
0.00 0.25 0.50 0.75 1.00
Ordered index
0.00 0.25 0.50 0.75 1.00
ordered index
Λ [~] λ
1.00
0.75
0.50
0.25
0.00
0.00 0.25 0.50 0.75 1.00
ordered index
Figure 20: The simulation result is the band in the picture above, and the curves represent Test a and b.
The top left is the QQplot of p-values, the top right is QQplot with both confidence bands from all tests, and the others are the QQplot of p-values with confidence band for each
test. 35 6 Discussion
The purpose of this study was to examine the relationship between global warming and extreme rainfalls in Taiwan. To address this question, Taiwan rainfall and temperature data from
TCCIP were analyzed. Specifically, the PGEV model was used to fit the annual maximum of daily rainfall to infer the patterns of extreme values and possible impact due to climite change in temperature. Results indicated that the intensity of extreme rainfall in Taiwan increases as the temper
ature rises, while the frequency remains constant. However, the effects of global warming on the frequency and intensity of extreme rainfalls varied by region. In the north and south Tai
wan regions, the frequency of extreme rainfalls changed with temperature, while in the central and east regions, the intensity of extreme rainfalls changed. The intensity of extreme rainfalls became higher as the temperature increased, according to the intensity plot. According to the
5% return value distribution, extreme rainfalls will happen more frequently in the future due to temperature continuously increasing. These findings have significant implications for water use and disaster prevention in Taiwan. Water in Taiwan mainly comes from the heavy rains brought by the mei-yu and typhoon season. If the temperature continues to rise in the future, Taiwan will face more frequent and intense extreme rainfall events, which could lead to serious floods and other types of disasters. In 2009, typhoon Morakot caused catastrophic damage in Taiwan and is recorded as the peak in the functional box plots. The study also goes through the entire process of analysis without 2009 data, and the result is similar to the analysis with 2009. For the post inference and visualization, Gaussian Process was used to spatially smooth the parameter estimates obtained from the fitted PGEV models and to construct the confi
dence bands for the p-value QQ plot. In particular, the latter task was implemented based on simulations using the Gaussian Process with Gaussian copula models. The confidence bands
36 show that all tests comparing with models with and without covariates are significant, and the
p-values from three tests are not significantly superior to others. 37