Risk Analysis of Gold Prices in Pakistan Using Extreme Value Theory

Extreme value theory (EVT) is useful for modeling the impact of crashes or situations of extreme stress on investor portfolios. EVT is mostly utilized in financial modeling, risk management, insurance, and hydrology.,e price of gold fluctuates considerably over time, and this introduces a risk on its own. ,e goal of this study is to analyze the risk of gold investment by applying the EVT to historical daily data for extreme daily losses and gains in the price of gold.We used daily gold prices in the Pakistan BullionMarket from August 1, 2011 to July 30, 2021.,is paper covers two methods such as BlockMaxima (BM) and Peak Over,reshold (POT) modeling. ,e risk measures which are adopted in this paper are Value at Risk (VaR) and Expected Shortfall (ES). ,e point and interval estimates of VaR and ES are obtained by fitting the Generalized Pareto (GPA) distribution. Moreover, in this paper, return-level forecasting is also included for the next 5 and 10 years by analyzing the Generalized Extreme Value (GEV) distribution.


Introduction
Gold is a familiar valuable metal for investment as compared to silver, platinum, and palladium. Gold has been widely used in many industries not only in the manufacturing of jewelry but also in financing and investing, medical and dentistry, manufacture of electronics, components of computers, and aerospace. In the financing and investing industry, gold is held in the form of gold bullion. Gold bullion is a valuable precious metal that is usually in two main forms: gold bullion bars and gold bullion coins. Nowadays, many institutions are including gold as a major asset of the world's investment due to their stability which plays an important role in the global economy. Gold is a good investment, neither in the shortterm nor the long-term, and it is a great asset that can be converted into paper money [1].
Gold is one of the oldest and most widely-held commodities used as a hedge against the risk of disruptions in financial markets, but gold prices fluctuated substantially over time: slowly and gradually, and this introduces the risk on its own. e gold prices in Pakistan are influenced by many factors, such as the value of the US dollar, gold production, consumer demand, devaluation of Pakistani currency, and inflation. ese factors which affect the gold prices are also discussed by Chaithep et al. [2].
In order to disclose the nature of the risks, under extreme situations, and finally, to avoid the risks to the most degree, we need certain risk measures. Extreme value theory (EVT) is for assessing the asymptotic probability distribution of extreme values for modeling the tail part of the distribution where the risk exists. EVT is playing an important role in dealing with modeling rare events. EVT applications are found in many fields where extreme events may occur, including hydrology [3], insurance [4], and finance [5]. EVT provides a solid framework for estimating the behavior of extreme events. e performance of EVT is better than the counterparts for predicting unexpected extreme movements because of focusing on the tails of the distribution [6].
Mainly, there are two broad methods of applying EVT; Block Maxima (BM) approach is based on the Generalized Extreme Value (GEV), while Peak Over reshold (POT) approach is based on the Generalized Pareto (GPA) distribution. e BM method is the primary and the oldest method of estimating extreme events in EVT. is method was given by Gumbal [7]. e BM method considers the distribution of the maximum order statistic in a GEV distribution that best fitted to the series of extremes observations. In this method, a single-highest value is taken from a block. e second highest value in this block may be larger than the highest of another block but this value is not accounted for. is is the conventional way of modeling extreme events and is used mostly in hydrology and other engineering applications. Further, the BM method is not suited for financial time series because of volatility clustering. Due to this phenomenon, extreme events tend to follow one another; the POT method is efficient as compare to the BM method because it uses the data more efficiently (it consider several large values instead of only the largest one).
is approach consists of using observations that are above a given threshold, called exceedances. erefore, the POT method has become the method of choice in financial applications. is method can be modeled by the Generalized Pareto (GPA) distribution [8].
In Figure 1(a), the observations X 2 , X 5 , X 7 , and X 11 show that the BM for four periods with three observations each. is is a traditional way of modeling extreme events and is used effectually in hydrology and other engineering applications. In Figure 1(b), the observations X 1 , X 2 , X 7 , X 8 , X 9 , and X 11 exceed the threshold (u), which represents extreme events. e distribution of exceedances above a threshold is based on a theory developed by Fisher and Tippett [9] that gives the limiting distribution of sample maxima provided the series which have certain mathematical properties [10]. e risk measures used in this paper are Value at Risk (VaR) and Expected Shortfall (ES). e VaR is a statistical technique that estimates the risk of an investment. In other words, VaR is the maximum possible loss or gain that could happen in investment over a given time period. ES is the average amount of all losses which are greater or equal to VaR. ere are two approaches, parametric and nonparametric, which are used to compute VaR and ES. In this paper, we used EVT that is a parametric approach for the estimation of VaR and ES.
Ren and Giles [11] applied the EVT in daily returns of crude oil prices in the spot market from years 1998 to 2006. ey concluded that the POT method by analyzing GPA distribution provides an effective risk measure such as VaR and ES. Jang [12] used the daily gold prices data from 1985 to 2006 to clear the main idea of EVT and discussed the tail behavior. He computed VaR and ES by using GPA distribution and normal distribution, and results show that the GPA distribution is superior to the normal distribution for VaR and ES estimations. Chaithep et al. [2] examined the EVT approaches to daily gold prices from the years 1985 to 2011 in the U.S. dollar. He concluded that VaR estimation by the GPA distribution model is more efficient than GEV distribution. Chen [13] analyzed the gold investment risk by using the EVT to historical daily data from years 1968 to 2014 in the London market. She compared the positive and negative return results of VaR and ES of different metals such as gold, silver, and platinum.
Chinhamu et al. [14] focused on the importance of the EVT approach for tail-related risk measures taking data from years 1969 to 2012. ey concluded that GPA distribution gives best results as compared to normal and Student's t distribution to estimate tail risk measures such as VaR and ES. In [15], the authors explained that the classical extreme value statistics consists of two fundamental techniques: the BM and the POT approach. It seems to be a general consensus among researchers in the field that the POT method makes use of extreme observations more efficiently than the BM method. POT is preferable to quintile estimation, while BM is preferable to return level estimation. In [16], the authors presented a study to assess the quality of VaR forecasts in various states of the economic situation. Two methods based on the EVT were compared, BM and POT. Forecasts were made on the daily closing prices of 10 major indices in European countries. Karahanoglu [17] analyzed the VaR method for one single portfolio or for all similar portfolios, and it hampers the opportunity for comparison in the finance market. e dataset is chosen as the period from January 2, 2014, to April 20, 2020, when all returns are recorded daily. e results show for similar portfolios, and different VaR methodologies and different back testing processes must be applied for the best fit. In [18], the author analyzed the VaR which has long been the standard risk measure in financial risk management. In this paper, he shows how EVT can be applied in estimating ES using an unconditional POTmethod. He measures the VaR as the standard measure of financial risk. In terms of ES, we find that cryptocurrencies are the high-risk asset and fixed income the shield. In [19], the authors explained EVT to model the tail behavior of the J580 Index for the right tail and the left tail. ey used the POTmethod to estimate the threshold using the Pareto Q-Q plot, Hills plot, and mean residual live plot. e MLE method was used to estimate the parameters. Extreme tail risk measures were estimated at relatively high quantile at the 95%, 99%, and 99.5% confidence intervals. e VAR and ES are used to examine the extent of the potential extreme downside and upside risk of the J580 Index. e results indicate the importance of the GPA distribution in fitting suitably the tails of financial time series data characterized by extreme events. e rest of the paper is organized as follows: In Section 2, we briefly describe the GPA and GEV distribution, Block Maxima, and Peaks Over reshold method for establishing extremes. Also, it describes the measures of extreme risks VaR and ES. Empirical results obtained in GPA and GEV distribution are discussed in Section 3. Finally, Section 4 concludes the study.

Basic Fundamental Assumptions.
Basically, it is compulsory to check the basic assumptions of data before conducting the analysis in modeling of any annual maximum series approach in the field of statistical hydrology. e basic and fundamental assumptions are stationary, randomness, homogeneity, and independence.

e Homogeneity Assumption.
e problem of homogeneity largely rises due to unexpected events in the observed sample data. is is occurring due to changes that occurred in the trend of the series. If a heterogeneity problem appears in the observed data series and frequency analysis results become unsuitable, then take out a particularly statistical investigation. In this situation, it is necessary to check the assumption of homogeneity before starting the statistical analysis of a given dataset. For checking the homogeneity assumption, we apply the Mann-Whitney U test. It is a nonparametric test that can also be used in the emplacement of an unpaired t-test and in the extreme type of data (such as gold price, rainfall, and flood data) that present variability that damaged the property of homogeneity.

2.1.2.
e Independence Assumption. e word independence means that the observed value of the sample is not the nonexistence of any other sample value. Puricelli [20] also discussed the historical hydrology data series such as the sequence of high and low, catchment flows, ARMS, annual maxima, or minima. e dependence is expected on the daily series, while no dependence is expected for periodic or annual series.
e independence assumption is formerly tested before performing the statistical analysis for hydrological variables. To check the independences' assumption, Wald and Wolfowitz [21] explained a test to check the independence assumption for the observed data. WW test is a nonparametric family test. e WW test is used to find whether the observed sample value is not affected by the existent of remaining sample data elements.

2.1.3.
e Assumption of Stationarity. e stationary assumption is important on sample data for their probability distribution, and related parameters are fixed with respect to time. In nonstationarity conditions, monotonic and nonmonotonic trends (unexpected changes at some point in time), we examined whether the sample data are stationary with respect to monotonic nonlinear or linear trends. ey help us in nonparametric tests based on Spearman's rank correlation coefficient.
e British psychologist Charles Spearman from 1863 to 1945 performed a nonparametric estimate of the statistical dependence between two random variables. It is shown by and named as in the Spearman rank correlation. ρ is statistical dependence, and it compresses a linear dependence as in usual Pearson's linear correlation coefficient.

e Assumption of Randomness.
NERC tests are used to check the randomness and find whether a dataset has an appreciable pattern and also whether the process that generated it is significantly random. If the sample of a hydrology variable can appear no randomness, then statistical dependence among its element has shown no homogeneity, and nonstationary causes can be also natural or man-made. Flood, climate fluctuations, and earthquakes are the natural causes that are mostly allied to where evolution causes are connected with land-use changes and building of big pool dams upstream, and humans make climate change.

Generalized Extreme Value Distribution.
e Generalized Extreme Value (GEV) distribution is three-parameter flexible distributions. It consists of Gumbal, Frechet, and Weibull distributions. e probability density function (PDF) of GEV distribution is given as where ξ is the shape parameter, μ is the location parameter, and α is the scale parameter. e shape parameter always shows the thinkness of the tail. So, if the shape parameter gives large value (means give positive value), then tail of the distribution is the thicker. Mathematical Problems in Engineering 2.3. Generalized Pareto Distribution. According to Pickands [22], the probability distributions of exceedances in an event are approximated to the Generalized Pareto (GPA) distribution. By introducing the shape parameter α and location parameter β, the two parameter GPA distribution has the following representation [23]. e PDF of GPA distribution is given below as Further, if the value of the shape parameter is high, then tail follows heavy-tailed distribution.

Excess Distribution.
For a random variable X, the excess distribution function F u above a certain threshold u is defined as where x represents the size of exceedances over threshold u. Further, if we denote F as the distribution function for X, then we may write A fundamental theorem in EVT, by Pickands [22], identifies the asymptotic behavior of these exceedances with Generalized Pareto (GPA) distribution. e excess distribution function F u can be well approximated by GPA distribution for large enough u: 2.5. Peak over reshold Method. We fit a Generalized Pareto (GPA) distribution to our dataset; we adopt the Peak Over reshold (POT) method that focuses on the distribution of exceedances above some high threshold. For x − u ≥ 0, we can rewrite the excess distribution function (4) as which allows us to apply the POT method: ere are two steps in applying the POT method. Firstly, we need to choose an appropriate threshold. Secondly, fit the GPA distribution function to data. Given the choice of a sufficiently high threshold, we may estimate F(u) by 1 − N u /n where n is the total sample size and N n is the amount of observations above the chosen threshold. And, F u (x − u) can be estimated by a GPA distribution using maximum likelihood estimation; we then can obtain the following tail estimator [11]: 2.6. Selection of reshold. e selection of the threshold level is an issue of adjusting bias and variance. In [24], the authors brought up that if a threshold is too low, it is more likely to violate the asymptotic property of the model and cause bias; if we take a threshold too high, it will generate few exceedances with which the model can be estimated and results in high variance. A basic procedure is to select a threshold as low as possible, as long as the limiting approximation of the model can give a reasonable result. So far, there is no algorithm-based method available for the selection of the threshold (u). Many researchers including McNeil [4], Danielsson and De Vries [25], and others have utilized this data-analytic issue, but none has provided a convincing solution.
In this paper, we used graphical methods for selection of the threshold level. In graphical methods, the most widely used plots are mean excess (ME) plot and shape parameter plot. ese plots are also used by Ren and Giles [11], Chen [13], Chinhamu et al. [14], Jakata and Chikobvu [15], and many other researchers for the selection of the appropriate threshold level. [26] introduced the simple mean excess (ME) plot. e first approach for threshold selection utilizes the ME plot. A mean excess function is the mean of the exceedances over a certain threshold (u). For a random variable X with right end point F x , the mean excess function is defined as

Simple Mean Excess Plot. Davidson and Smith in
For u < F x , if the underlying distribution of X > u has a GPA distribution, then the corresponding mean excess is In equation (8), we can clearly see that the mean excess function must be linear in u. More precisely, it follows a GPA distribution if and only if the mean excess function is linear in u [24]. is gives us a way of selecting an appropriate threshold.
Given the data, we define the empirical mean excess function as where n is the sample size. e empirical excess plot is a graphical representation of the locus of (u, e n (u)), and we can examine this plot to choose the threshold u such that e n (u) is approximately linear for X > u. We use R software to make the plot. More discussion will be presented in Section 2.6.

4
Mathematical Problems in Engineering

Parameter Estimation.
ere are several techniques for estimating the parameters such as maximum likelihood estimation (MLE), method of moments, and method of probability-weighted moments.
In this paper, the MLE method is used for estimating the parameter of GEV and GPA distribution because the MLE is asymptotically normal and allows simple approximations for standard errors and confidence intervals.
Suppose that block maxima X 1 , X 2 , X 3 , . . . , X k are independent from a GEV distribution, and the log-likelihood function for the GEV distribution from the equation (1) is given as Smith [27] declared that, for ξ > 0.5, the MLE for ξ, u, and σ satisfies the regular condition and therefore having asymptotic and consistent properties. e number of blocks k and the block size form a crucial trade-off between variance and bias of parameter estimation [2].
Given that we have a adequately high threshold u and assuming there are m observations with then, the logarithm of the PDF of x i can be derived from equation (2) as Hence, the log-likelihood function L(α , β | x i − u) for the GPA distribution is the logarithm of the joint density of the m observations is as We can obtain the estimates for α and β by maximizing the log-likelihood function of the subsample under a suitable threshold (u) [14].

Assumption of Poisson Distribution.
e POT approach has many features in the modeling of exceedances in which the observations exceed or occur over a given period of time.
In the POT method, the general and simple rule is that the numbers of exceedances follow a Poisson process which is however not necessary for modeling the BM approach [28]. However, Rosbjerg et al. [29] considered this assumption compulsory for modeling of the POT method. A Poisson distribution is given in the equation below: where e is the base of natural log, μ is the number of success mean, and x is the successes number.
In this assumption, the occurrences must be independent and the probability of occurrence of an extreme event at a time "t + s" does not depend on another extreme event that occurs at a time "t." is assumption may not fulfill in the case of natural processes (i.e., outbreak of a disease, floods, hurricanes, and draught) as they occur in clusters. In this condition, the decluttering method is normally used. e property of Poisson distribution is that it has same mean and variance. However, the gold price dataset usually contain some outliers, which may result in violating the Poisson process assumption. Cunnane [28] suggested a criterion on the bases of confidence interval to test this assumption; the test is called as "Fisher Dispersion Index" (FDI). e test statistics proposed by Cunnane is based on approximating the Poisson variants through normal variety. e number of peaks that occur in the i th year is expressed as n i and is normally distributed with mean λ and standard deviation λ; then, the test statistics: where n i is the number of exceedances that happen in the i th year and λ is the mean of the numbers that exceed the threshold level. Following the chi-square distribution with degree of freedom b � N-1, where N is total number of years, this is valid for N greater than 5 (N > 5). e test statistics that perform at "α" significance level under the null hypothesis H 0 is accepted if the value of R lies between  [30].
For a random variable, X is usually the return in some risky financial market with distribution function F over a specified period, and the VaR for a given probability p can be defined as the pth quartiles of F is as where F − 1 is the quartile function. VaR is a common measure of extreme risks, and we used GPA distribution to approximate this measure. In particular, using (6), we obtain VaR as where β and α are the maximum likelihood estimates of the GPA distribution parameter.

Expected Shortfall. Expected Shortfall (ES) is proposed
by Artzner et al. [31], and it is the average amount of loss, or in other words, ES is the average of all losses which are equal to or greater than the VaR level.
In contrast to VaR, ES measures the risk of the market by considering both the size and likelihood of the losses above a particle threshold. ES gives the expected size of return that exceeds VaR for a probability level p: and equivalently, where the second term above represents the excess distribution F VaR p (x), treating VaR p as the threshold level, proceeding as before. If the threshold VaR p is sufficiently large, then F VaR p (x) is GPA distribution i.e., us, the mean of the excess distribution F VaR p (x) can be calculated as where α < 1 and substring into equation (19) will yield

Return Level Estimation.
e most widely useful application of the block maxima approaches is estimating the return levels based on the Generalized Extreme Value (GEV) distribution. It is the quantile points of the GEV distribution: where M is the distribution of maxima observation of equal length over a successive nonoverlapping period. e expected level to be exceeded in one out of k periods of length n where k is the period and n is the length of this period. e return level can be used as a measure of the maximum loss of a portfolio, a rather more conservative measure than the Value at Risk [32].

Data Description.
e data used in this paper are the daily gold prices in the Pakistan Bullion market over the time period of August 1, 2011 to July 30, 2021, and the data is taken from the following website: https://www.bullion-rates. com/gold/PKR/2007-2-history.htm. e daily gold returns are calculated by taking log differences of the daily gold prices: e total number of observations is 2707, and there are 2706 gold returns, among which 1304 are negative returns and 1402 are positive returns. We separated the positive returns and the negative returns and then take the absolute values of the negative returns because Extreme Value eory (EVT) is only defined for nonnegative values [14]. Table 1 provides a summary of descriptive statistics for the considered return series.
In Table 1, we observed that the mean of daily return is 0.00012 which is a positive indication that the overall gold prices are increasing during the considered time period. e standard deviation (SD) is 0.0128, the maximum value is 0.0871, and the minimum value is − 0.1172 in daily return. e magnitude of the average return is very small as compared to the standard deviation. Further, the large kurtosis of 10.3284 indicates the distribution of returns has a fat tail.
e Jarque-Bera statistic has zero probability; it indicates that we can reject the null hypothesis in favor of the alternative hypothesis, so the distribution is nonnormal [13]. e daily gold prices in the Pakistani rupee are shown in Figure 2, and daily returns are shown in Figure 3. e plot of the daily gold prices (Figure 2) shows a substantial increase since 2011 with lots of fluctuations, and the graph of daily returns ( Figure 3) shows a large number of extreme values in the dataset. According to [12], the graph of daily return confirms the volatility of the Pakistan Bullion Market. e descriptive statistic and time series plots are obtained with the help of Eview 8 software.
e autocorrelations (ACF) pattern helps us to show the unpredictability in returns and volatility clustering. Autocorrelations investigate the correlations between today's value and value in the past time series. Figure 4 gives ACF of return series, and Figure 5 gives the ACF plot of squared returns series. e return autocorrelations are almost all insignificant while the squared returns have significant autocorrelations. Furthermore, the squared return autocorrelations are all positive which is highly unlikely to occur by chance. e figures give significant evidence for both the unpredictability of returns and volatility clustering [14].
Further, to check the stationary of the daily return of gold prices, we fit Augmented Dickey-Fuller (ADF) and Philips Peron (PP) Unit Root tests' results, which are given in Table 2.
In Table 2, we are applying the ADF and the P-P unit root tests on the return series of the gold prices to show the stationarity. In the ADF test, we used Schwartz Information Criterion to set at lag 0, and Bartlett Kernel spectral estimation method was used in the P-P test. e null hypothesis of unit root is rejecting in both tests (p ≤ 0.001); in this situation, we accept the alternative hypothesis and conclude that the daily return of gold prices series are stationary.

Basic Fundamental Assumptions.
It is compulsory to check the basic assumptions of data before conducting the analysis; the basic and fundamental assumptions are stationarity, randomness, homogeneity, and independence. e details of these assumptions are available in Section 2.1.
In Table 3, the results are satisfying the basic and fundamental assumptions of stationarity, randomness, homogeneity, and independence. We checked all fundamental assumptions for monthly, quarterly, and yearly Block Maxima (BM) of negative return. e test statistic value for each BM of a negative return is small and p ≥ 0.050. erefore, we concluded that the monthly, quarterly, and yearly BM of a negative return is satisfying the fundamental assumptions of stationarity, randomness, homogeneity, and independence.

Determination of Best-Fitted Distribution.
In this paper, first, we used the Block Maxima (BM) approach in which the Generalized Extreme Value (GEV) distribution is fitted. e same distribution was used by Pratiwi et al. [33]. In order to find the best fitting of GEV distribution, both graphical and formal tests of goodness-of-fit are applied. In graphical, to find the best-fitted distribution for BM of negative return, we used crude residual plot and Quantile-Quantile (Q-Q) plot. In Figure 6(a), crude residual and Q-Q plot are for monthly BM of negative return, Figure 6(b) for quarterly BM of negative return, and Figure 6(c) for yearly BM of negative return. e crude residual plot for monthly, quarterly, and yearly BM of negative return shows the concave pattern from the straight line, and it indicates that the data series come to a heavy-tailed distribution. e Q-Q plot looks linear to show that the GEV distribution is best-fitted for BM of negative return.
After graphical view, we concluded that GEV distribution is best fitted for monthly, quarterly, and yearly BM of negative return. Now, Kolmogorov-Smirnov (K-S), Anderson-Darling (A-D), and chi-square tests are applied to find the best-fitted distribution for BM of negative returns.
In Table 4, some formal tests K-S, A-D, and Chi-Square tests' goodness-of-fit test are applied to find the best-fitted distribution. e test-statistics value of these goodness-of-fit tests is calculated by the R software. Table 4 shows the teststatistic value of GEV distributions and their p value. According to the K-S, A-D, and the Chi-Square test, all the p values are greater than 0.05, so GEV distribution is the bestfitted distribution for monthly, quarterly, and yearly BM of negative return.
We are focusing on GEV distribution which is a bestfitted distribution for BM of negative return. Next, we calculate the parameters of the selected distribution.

Parametric Maximum Likelihood Estimates.
e Block Maxima (BM) of negative returns has been fitted to the Generalized Extreme Value (GEV) distribution with three different block sizes (yearly, quarterly, and monthly). We estimate the parameter by using the Maximum Likelihood (MLE) method. Refer to Section 2.7, equation (11) for the maximum likelihood function. e R software is used for modeling the extreme events. Table 5 shows the the estimates of three-parameter location, scale, and shape with their standard error of GEV distribution for three different blocks size (yearly, quarterly, and monthly). According to Jang [12], all block numbers show   Mathematical Problems in Engineering that the shape parameter is positive. All the parameters have decreasing estimated standard errors as the number of blocks increased. In general, the BM of the negative returns follows a Frechet family of GEV distribution for yearly, quarterly, and monthly block frames. e Frechet type of GEV distribution confirms that the original series has a fat tail.

Estimation of Return Level.
In this paper, we divide the n sample block into monthly, quarterly, and yearly blocks by taking the highest value in each block. e point estimate is estimated for the sample data (obtained through the BM approach) based on the MLE method using the probabilities of nonexceedances for different return periods. By putting the nonexceedances probability to the distribution's, the certain return period is computed. Refer to Section 2.8 for the maximum likelihood function. e R software was used for modeling the extreme values. In Table 6, we estimate the point and interval estimates for the next five years and ten-year return period. We have monthly BM of the negative returns, and for the next five years' returns period, the point estimate is 0.0293, the corresponding probability of nonexceedances is 0.80, and the probability of exceedances is 0.20. It means that there are 80% chances that monthly maximum loss will not exceed 2.93% once in the next five years on average and a 20% chance it will be exceeded.
For the next ten-year return period, the point estimate is 0.0363 and the corresponding probability of nonexceedances is 0.90, and the probability of exceedances is 0.10. It means that there are 90% chances that the monthly maximum loss will not exceed the value of 3.63% once on average in the next ten years and a 10% chance it will exceed. According to Jang [12], the point estimate is increased with decreasing block sizes. With the increasing time interval, the point estimate increases, and the confidence interval is wider. It reveals that time is a risk factor.
Similarly, we interpret the return level for quarterly and yearly BM of the negative returns.
3.6. Selection of reshold. As we have discussed in Section 2.6, there is no algorithm-based method available for the selection of the threshold; therefore, we used a graphical method for the selection of the threshold level. Two plots could help with the determination of the threshold level. Figure 7 presents the Mean Excess (ME) plot, and Figure 8 presents the estimate of the shape parameter for positive returns, and Figure 9 presents the estimate of the shape parameter for negative returns. ese plots are obtained with the help of R software. e ME plot is helpful in detecting graphically the quintile above which Pareto's relationship is valid. In Section 2.5, the ME plots are approximately linear in the threshold (u) given that the underlying distribution of sample data is a Generalized Pareto (GPA) distribution. More specifically, the ME plot of the data can be used to distinguish between light and heavy-tailed models; the plot of a heavy-tailed distribution shows an upward trend, a medium tail shows a horizontal line, and the plot is downward-sloped for lighttailed data. Common ground in our sample data is that both the ME plots of positive and negative returns have an upward trend followed by an irregular the portion in the far end. Figure 7 shows ME plots for positive and negative returns. According to Chen [13], we take the threshold (u), where the ME plots are linear and upward-sloping. In this paper, the ME plots' upward-sloping is started from 0.001 to 0.015 thresholds' level. erefore, we take the range of threshold level where the plots look like linear and upward slope in different segments. Figures 8 and 9 show that the shape parameter plots of MLE estimates the shape and scale parameters under different thresholds' level. e 95% level of confidence intervals is computed for upper and lower dashed lines of plots. According to Chen [13], we took the threshold level before which the shape parameters are stable. According to Ren and Giles [11], we took the threshold level before the plot is comparatively flat. e gain shows an approximate upward trend in the threshold from 0.001 to 0.015. e loss shows approximate linearity with a slightly upward trend in the threshold from 0.001 to 0.010. erefore, there is some evidence to choose thresholds from 0.001 to 0.012 for the right tail and from 0.001 to 0.009 for the left tail.

Parameter Estimation.
Given the thresholds selected in the previous section, we could estimate the shape and scale parameters of the corresponding GPA distribution. Refer to Section 2.7, equation (13) for the maximum likelihood     function. e R software was used for modeling extreme events. Table 1 summarizes all the parameter estimates and their standard errors under different thresholds.
In Table 7, the confirmed choice of the threshold is u � 0.008 for positive returns and u � 0.001 for negative returns. For positive returns, there are 99 observations above the threshold which is 0.008. Compared to the sample size 1402, the tail only accounts for 7.06% of the total distribution. As mentioned in previous sections, EVT only analyzes the extreme events; thus, even if we started with a large sample, we contend with few observations. Moreover, compared to using weekly returns, daily returns provide a large sample size. is reduces the variance by increasing the number of observations in the tail. It is the same situation for negative returns.
Since the maximum likelihood estimator is asymptotically normal, the associated MLE of the parameters under selected thresholds is statistically significant at the 5% significance level [13].

Assumptions of Poisson Distribution. According to
Cunnane [28], if the dataset contains outliers, the assumption is tested through Poisson distribution statistics with fixed confidence intervals. He also shows the use of the criterion of the Fisher Dispersion Index to check the Poisson   process assumption. e details are discussed about Poisson assumption in Section 2.8. e calculation is done in Excel by using equation (15). In Table 8, the test statistic R is tested at 1% and 5% significance levels. For positive return, the threshold 0.008 and 0.012 of a sample of 99 and 39, respectively, and observations are considered for positive return. e test statistic R is 15.44 and 11.51 which falls in the acceptance region for both 1% and 5% significance levels which indicates that the sample obtained through the POT approach follows the Poisson process assumption. e sample of sizes 99 and 39 is selected on the level of threshold. is sample satisfies the Poisson process assumption; therefore, further POT modeling analysis for positive return analysis is carried on the basis of this sample.
For negative return, the threshold is 0.001 and 0.003 of a sample of 962 and 493, respectively, and observations are considered for a negative return. e test statistic R is 10.14 and 17.89 which falls in the acceptance region for both 1% and 5% significance levels which indicates that the sample obtained through the POT approach follows the Poisson process assumption. A sample of sizes 962 and 493 are   selected on the level of threshold. e samples satisfy the Poisson process assumption for a negative return. erefore, further POT modeling for negative return analysis is carried on the basis of these two sample data.

Assumption of Generalized Pareto
Distribution. e exceedances in the POT approach are considered to identically distribute with exponential distribution, but in the presence of outliers, the exponential distribution fails to fulfill the model requirements. In recent years, GPA distribution attained high consideration in POT modeling because of its flexibility. We test the goodness-of-fit of GPA distribution above the threshold level by using K-S and A-D tests. e results of GPA distribution of goodness-of-fit for positive and negative threshold are given in Table 9. Table 9 shows goodness-of-fit results for GPA distribution of positive and negative daily return at the threshold level; for positive return, the value of test statistic at threshold 0.008 is 0.1105 and 0.3312, respectively, and the p value are 0.1653 and 0.9128 for both K-S and A-D tests. Since p ≥ 0.050, therefore we can say that the exceedances or the POTseries follow the GPA distribution. e similar behavior is observed for the other threshold level.
It means that daily positive return of gold prices at threshold level 0.008 and 0.012 and daily negative return of gold prices at threshold level are 0.001 and 0.003 which are best fitted for GPA distribution.
EVT suggested that the excess distribution above a suitable threshold of daily returns should follow a GPA distribution. To determine how the GPA distribution fits the tails of the return distribution, we plot the empirical distribution of exceedances along with the cumulative distribution simulated by a GPA distribution and compare the results visually. Figures 10-13 are the residual diagnostic check of GPA distribution fit to daily positive and negative returns of gold prices. Figures 10 and 11 provide the plots of excess distribution plot and tail of underlying distribution plot for gains and loss, and the plots show that the GPA distribution fitting to 99 exceedances at the threshold u � 0.008 with the shape parameter estimate 0.1282, and the GPA distribution fitting to 39 exceedances at the threshold u � 0.012 with the shape parameter estimate 0.0978. For losses, the plots give the GPA distribution fitting is 962 exceedances at the threshold u � 0.001 with the shape parameter estimate of 0.1251, and the GPA distribution fitting is 493 exceedances at the threshold u � 0.003 with the shape parameter estimate 0.1762. For both positive and negative returns, the graph of the excess distribution function follows closely the trace of a corresponding GPA distribution, implying that the GPA distribution models the exceedances very well [11]. Figures 12 and 13 provide the plots of scatter plot of residuals and Q-Q plot for gains and loss and are the quantile plots to measure the quality of best-fitted GPA distribution. e GPA distribution is sensibly fitted to the exceedances above the threshold when the Q-Q plot shows a linear pattern. If the scatter plot of residuals does not show any pattern, it means that the peaks values are independent [14].
Overall, based on different fittings with different values of the threshold and associated parameter estimates, the GPA distribution models the tail behavior of our daily returns very well and the fits exhibit reasonable robustness to the choice of thresholds.

Risk Measures' Estimation.
e Value at Risk (VaR), Expected Shortfall (ES), and their respective confidence intervals for both positive and negative returns give different levels of the threshold and different tail quintiles. e risk measures used in this paper are VaR and ES. e VAR is a statistical technique that estimates maxima possible risk over one trade day and ES is the average amount of all losses which are greater or equal to VaR. ere are two approaches, parametric and nonparametric, to compute VaR and ES. In summary statistics, we suggest that   the returns' distribution has a fat tail. In fat tail distribution, the nonparametric approaches for computing VaR and ES failed. erefore, in this paper, we used EVTfor estimation of VaR and ES which is a parametric approach [14].
In this paper, Value at Risk (VaR) measures the best and worst-case scenario on the Pakistan Bullion Market value of the par-gram gold price over one trade day, given a specified degree of confidence. We first consider the cases of point estimates under the lower threshold for both tails (0.008 for the right tail and 0.001 for the left tail) with statistics shown in Table 10. e calculated VaR is 0.0325 with 1% probability (99th percentile) for the right tail. at is, given usual conditions, we expect a daily gain in the value of par-gram gold prices in the Pakistan Bullion market would not increase by more than 3.25%, and the average gain above this level will be 4.06%. In other words, the market value, with a probability of 1%, would be expected to gain by Rs.32,500 or more, and on the average basis, the gain is by Rs. 40,600 if we have an investment of 1 million Pakistani rupees in that market.
On the contrary, VaR is estimated as 0.0185 with a probability of 1% (99th percentile) for the left tail. e worst daily loss in the market value could fall below 1.85%, and the average loss below this level will be 2.42%. Put differently, if we invest 1 million Pakistani rupees in the Pakistan bullion gold market, we are 99% confident that our daily loss at worst will not exceed Rs.18,500, and on the average basis, the loss is Rs.24,200 during one trade day.
On the contrary, VaR is estimated as 0.0185 with a probability of 1% (99th percentile) for left tail. e worst daily loss in the market value could fall below 1.85%, and the average loss below this level will be 2.42%. Put differently, if we invest 1 million Pakistani rupee in Pakistan bullion gold market, we are 99% confident that our daily loss at worst will not exceed Rs.18,500, and on the average bases, the loss is Rs.24,200 during one trade day.
Similarly, at a lower quantile of 95-level, the estimated VaR is 0.0223 for gains and 0.0112 for losses. It can be stated at, with the 95% confidence, the expected market value of par-gram gold prices would not gain by more than 2.23% for the best-case scenario or lose more than 1.12% for the worstcase scenario within the one-day duration. e adequacy of VaR estimates' essential depends on their accuracy. e easiest way to examine this method is to construct the confidence intervals. For instance, we discuss the right tail here. An approximate 95% confidence interval constructed for the 99-level quintile VaR is (0.0257, 0.0541) under threshold is 0.008 and (0.0289, 0.0579) under threshold is 0.012; for the 95-level quantile VaR, the associated 95% confidence intervals are (0.0193, 0.0284) and (0.0232, 0.0425).  Some characteristics of the estimation can be summarized as follows: (1) under different thresholds, the estimates of VaR exhibit strong stability; (2) given the quantile level, under either the lower or higher threshold, the corresponding VaR estimate in the right tail is larger than in the left tail, but the difference is small, implying that the behavior difference in both tails is likely to be small.

Conclusion
e high unpredictability of gold prices in Pakistan bullion markets requires the application of effective risk management. Extreme Value eory (EVT) is a capable method to estimate the effects of extreme events in risky markets. In this paper, we have described the EVT which can be used to model the tail-related risk measures, such as Value at Risk (VaR), Expected Shortfall (ES), and return level by applying it to the daily returns of gold prices in the Pakistan Bullion Market. We used the POT method to estimate the threshold using the mean excess plot and shape parameter plot. e MLE method was used to estimate the parameters of GEV and GPA distribution. e point and interval estimates of VaR and ES calculated under different high quantile levels show strong stability through a range of the selected thresholds, suggesting the accuracy and reliability of the estimated quantile-based risk measures.
Moreover, it is included in the paper return level forecasting in the next 5 and 10-years by analyzing the Generalized Extreme Value (GEV) distribution. It is to be noted that, with the increase in the time interval, the point estimate increases slowly and the confidence interval becomes wider. erefore, this indicates that time is a risk factor. Our application studies the heavy-tailed behavior in daily returns and the asymmetric characteristics in distributions, suggesting we treat positive besides negative returns separately. e results show that the EVT is a powerful technique to estimate the effects of extreme events in financial markets.
Our EVT-based VaR approach provides quantitative information for analyzing the point of potential extreme risks in bullion gold markets. Organizations and corporations may use this technique as a means of risk management. Further, for the people who want to invest in the Pakistan bullion markets, the estimates of VaR and ES provide quantitative indicators for their investment decisions.

Data Availability
We used the daily gold prices in Pakistani rupees from Bullion Rates > Gold Price History in Pakistan Rupees (PKR) for February 2007 from the following website: https:// www.bullion-rates.com/gold/PKR/2007-2-history.htm.

Conflicts of Interest
e authors declare that they have no conflicts of interest.