Markovian approach to the frequency of tropical cyclones and subsequent development of univariate prediction model

A tropical cyclone is one of the most devastating meteorological events. In recent years, we faced some very severe cyclones to a super cyclone successively that caused heavy damages to life and property during the helpless situations of the global pandemic. In this paper, we studied the frequency of cyclones from the years 1891 to 2019, i.e., for 129 years, on the Arabian Sea Basin, Bay of Bengal Basin, and land. We have categorized the cyclones according to their wind speeds: (i) cyclonic storms and severe cyclonic storms (CS + SCS) and (ii) depressions, cyclonic storms, and severe cyclonic storms (D + CS + SCS) where depressions, cyclonic storms, and severe cyclonic storms have wind speeds of more than or equal to 17 knots, 34 knots, and 48 knots respectively. We examined the Markovian dependence of the discretized time series of the two categories mentioned earlier for the first, second, third, and fourth order of a two-state Markov chain model. It is found that CS + SCS represents the first-order two-state (FOTS) model of Markov chain and D + CS + SCS represents the second-order two-state (SOTS) model of Markov chain. Thereafter, we have developed autoregressive models for the two categories and checked their goodness of fit using Willmott’s indices of orders 1 and 2. It is found that CS + SCS best represents the autoregressive model of order 5, whereas D + CS + SCS could not be efficiently represented by the developed autoregressive models. So we further developed autoregressive neural networks for D + CS + SCS and obtained a significant hike in the prediction yield. Nevertheless, it is found that both categories are not serially independent.


Introduction
Tropical cyclones (TCs) are synoptic-scale phenomena where a large mass of air swirls around a low pressure, counter-clockwise direction in the Northern Hemisphere and clockwise in the Southern Hemisphere. As referred to as the heat engine, a tropical cyclone is mainly fuelled by the latent heat of the moist air rising from the ocean. Gray (1968Gray ( , 1979 identified the six parameters necessary for the formation of tropical cyclones which are widely accepted even today. It consists of cyclonic low-level relative vorticity, a large value of relative humidity in mid-troposphere, conditional instability through a deep tropospheric layer, warm and deep oceanic mixed layer, weak tropospheric vertical shear of the horizontal wind, and location of disturbance a few degrees poleward of the equator. Varotsos et al. (2019) designed a cluster algorithm to predict the instability of the atmosphere as it is an important precursor for the formation of tropical cyclones. The information-monitoring-tracker proposed in this paper gave promising results.
Around 80 TCs are formed globally every year (Emanuel, 2003). The North Indian Ocean (NIO) contributes about 7% of the global TCs (Gray, 1979). These cyclones are primarily originating in the Bay of Bengal (BOB) and the Arabian Sea (AS). The stages of cyclones are classified as depression (D), deep depression, cyclonic storms (CS), severe cyclonic storms (SCS), very severe cyclonic storms, extremely severe cyclonic storms, and super cyclones based on their associated wind speeds. On average, there are about 5 cyclones every year in the NIO. For every 4 cyclones formed in the Bay of Bengal, 1 cyclone is formed in the Arabian Sea (Mohapatra et al., 2017). In recent years, we see an increase in the intensity of pre-monsoon cyclones in the AS. Rajeevan et al. (2013) suggested that epochal variation in the intensity of TCs over the AS is correlated with the epochal variation of vertical wind shear and TC heat potential. Mohanty et al. (2012) further pointed out that after 1950, the frequency of severe cyclonic storms has increased significantly by 71% in the BOB and 300% in the AS during the post-monsoon months. Deshpande et al. (2021) showed that there is a 52% increase in the frequency of cyclonic storms and a 150% increase in the frequency of very severe cyclonic storms for the epoch 2001-2019 in the Arabian Sea Basin. However, there is a decrease in cyclonic storm frequency by 8% in the Bay of Bengal for the same time period. Cracknell and Varotsos (2011) questioned whether the climate models should be deterministic, chaotic, or stochastic (Chattopadhyay and Chattopadhyay 2008). Mooley (1980) indicated that the intensification of storms into severe storms is distributed by binomial probability model, whereas the formation and landfall is based on Poisson stochastic processes. Mooley (1981) and Mooley and Mohile (1984) showed that the mean annual frequency of severe storms has increased over the NIO in the period 1965-1980at a 10% level. Srivastava Sinha De (2000 reported that in the period 1891-1997, there is a significant decreasing trend in the number of cyclonic storms in the NIO at 99% level which might be due to the weakening of Hadley circulation. This decreasing trend was more in the BOB than in the AS. Although the annual frequency was recorded as having a decreasing trend of nearly 15% per hundred years, there was an increasing trend of cyclones in November and May mainly contributed by the BOB (O.P. Singh et al., 2001). This paper also stated that during November there was an increase of 20% per hundred years in the rate intensification of cyclonic disturbances to the severe cyclonic stage. Later, using the regional climate model HadRM2, O.P Singh (2007) indicated that the climate change due to an increase in the atmospheric greenhouse gas concentration is the cause of this increasing trend during intense cyclone months May, October, and November. His simulation experiments also showed that the frequency of post-monsoon tropical disturbances in the Bay of Bengal will increase by 50% by the year 2050. It is a hot topic of research whether the increase in cyclonic frequency is related to an increase in SST. Even though Gray(1979) and Sikka (1977) showed that sea surface temperature (SST) of more than 26 °C is a parameter for the genesis of a cyclone, Pattanaik (2005) and Chan (2007) reported that the variability in the planetary-scale atmospheric circulation is the main cause of the interdecadal variability of cyclonic activity over the Indian region rather than the variability of SSTs over the region. Goldenberg et al. (2001) related the sudden increase in hurricane activity during the period 1995-2000 with the increase in SST and decrease in vertical wind shear in the North Atlantic. They also claimed that this increased activity will continue for the next 10 to 40 years. Later, Webster et al. (2005) also showed that the number of tropical storms and hurricanes cannot be correlated with increasing SST. However, it was shown in a single storm simulation that SST and its gradient played an important role in the peak intensity and track of a tropical cyclone (Mandal et al., 2007). Over the North Atlantic Ocean, it is found that the frequency of hurricane count has increased significantly in the last century  making two abrupt shifts in 1925-1926 and 1987-1988 which coincides with the increase in SST in that time interval (Varotsos et al., 2015). A strong seasonal link has also been shown between SST and El Niño/Southern Oscillation (ENSO) in the northern hemisphere, the SST time series being stochastic (Varotsos, 2013, Varotsos, 2014. Nina Črnivec et al. (2016) showed that the intensity of a cyclone also depends on the latitude when SST is changed. It states that intensification is more dependent on SST for higher latitude (say 25°N) than a lower latitude (say 10°N). Mandke and Bhide (2003) pointed out that during 1958-1988, the frequency of cyclones over the BOB decreased even though the SST increased. Nolan and Rappin (2008) also stated that on increase in SST in a radiative-convective environment, the wind shear prevents cyclone genesis. Gao and Zhai (2016) suggest that the latent heat flux and sensible heat flux play an important role in the intensification of the tropical cyclone over the western North Pacific Ocean. Hallam et al. (2021) showed that the intensity of the cyclone is most closely correlated with the ocean temperature profile of the upper 80 m of the ocean water and not the upper 50 m which was thought earlier. Lengaigne et al. (2018) showed that the air-sea interaction acts as negative feedback dampening the intensity of tropical cyclones over the Indian Ocean due to the induced cooling.
The frequency of cyclones is also speculated to be affected by various other phenomena or parameters. Gray (1984) correlated El Niño/Southern Oscillation (ENSO) with tropical cyclone activity in the North Atlantic, whereas Chan (1984) observed the same in the Northwest Pacific.
The influence of the Madden-Julian Oscillation (MJO) over the Australian region was studied by Hall et al. (2001) who showed increased cyclonic activity in the active phase of MJO which was strengthened during El Niño. Further Klotzbach (2014) also showed that TC activity is enhanced during and immediately following the active convective phase of the MJO while it is suppressed during and immediately following the convectively suppressed phase throughout the globe. B Kumar et al. (2011) related the decreasing trend of CS and SCS in the pre-monsoon with the increasing SST over the NIO in general and the BOB in particular, whereas in the post-monsoon season the frequency of tropical systems is positively related with Southern Oscillation Indices (SOI) and inversely correlated with MJO index for the period 1891-2008. This paper also mentions that there is an influence of El Niño and La Niña on the frequency of tropical cyclones over the BOB. Ng and Chan (2012) point out that the tropical cyclone activity is less in the El Niño year and more in the La Niña year during Oct-Dec and the possibility of ENSO and Indian Ocean Dipole (IOD) influencing tropical cyclone genesis and development due to variation in the atmospheric dynamic and thermodynamic conditions during the post-monsoon season. Deo and Ganer (2014) also explained the increase in cyclones over the AS due to variability in SST, wind shear, and energy metrics like Accumulated Cyclone Energy and Power Dissipation Index. They also showed that the cyclone season length, i.e. the number of days from the start of the cyclone season to the end of the cyclone season, is increasing at the rate of 0.33 days per year mainly due to pre-monsoon season length. Later, Baburaj et al. (2020) showed that there is epochal variability of cyclone frequency in the AS, whereas in the BOB there is a decrease in cyclone frequency in all three epochs which is related to epochal decadal variability in the equatorial Indian Ocean SST and vertical variation of the thermal profiles during the three epochs due to the warming of both atmosphere and ocean. However, other factors influencing the frequency of cyclones in the NIO might also be present which are yet to be discovered. Ki-Seon Choi et al. (2009) developed a multiple linear regression model that showed that the Tibetan Plateau snow cover is an important influence in the formation of tropical cyclones in Korea. Gray (1975) introduced the new factor Yearly Genesis Parameter to measure tropical cyclogenesis which was further modified by Royer et al. (1998) by including the impact of greenhouse gases.
The occurrence of cyclones causes havoc in the environment. Emanuel (2005) defined the term Power Dissipating Index based on the power dissipated by a hurricane in its lifetime that causes destruction, and his study for the period 1970-2004 shows that increase in this index over North Atlantic plus western North Pacific is partly due to increase in SST. Lovejoy and Varotsos (2016) showed that the response of the climate system has had a linear relationship with the radiative forcing for the last 50 years. The highresolution climate model analysis by Gupta and Jain (2020) suggests that in the twenty-first century, the frequency of the most intense cyclones in the BOB and AS will likely increase due to warming while the total number of cyclonic disturbances should decrease. Using atmospheric general circulation models under the Intergovernmental Panel on Climate Change (IPCC) A1B scenario and phase 5 of the Coupled Model Intercomparison Project (CMIP5) models under the representative concentration pathway (RCP) 4.5 and 8.5 scenarios, Murukami et al. (2014) indicated a decrease in the projected frequency of cyclones in the basins of the Southern Hemisphere, Bay of Bengal, western North Pacific Ocean, eastern North Pacific, and the Caribbean Sea and increases in the Arabian Sea and the subtropical central Pacific Ocean. Earlier, Emanuel (2013) showed that, on downscaling tropical cyclones of CMIP5 models for the period 1950-2005 and comparing with the twenty-first century, there is a noticeable increase in cyclone activity as well as an increase in the intensity of cyclones in the North Pacific, North Atlantic, and South Indian Oceans.
In recent times, the track of the cyclone can be comprehended 48-72 h in advance. However, the intensity and accompanying storm surge need more advanced models to be predicted more precisely (Sikka, 2006). S. K. Dube et al. (2009) showed the developments in storm surge prediction in the BOB and AS. Krapivin et al. (2012) designed a procedure using the percolation approach to forecast the beginning of a tropical cyclone earlier than traditional methods when the ocean-atmosphere system shifts from a stable background state to a state with a high level of instability. Being a major threat to life and property, its frequency, intensity, and track need to be properly estimated such that the required precautions can be taken and the budget for relief funds can be decided. Hence, it is a major concern for the IMD to analyze the frequency as well as the intensity of cyclones over the North Indian Ocean. Efstathiou and Varotsos (2012) showed that there is a persistent long-range correlation of the Sahel precipitation anomalies for the period 1900-2010 for all time lags between 4 months and 28 years. Later, Varotsos and Efstathiou (2013) have analyzed the long-term memory effect of the tropical cyclones using detrended fluctuation analysis and concluded the absence of any long-term correlation in the tropical cyclone frequency time series for the period 1870-2006 in the Atlantic. This paper also reveals the presence of white noise in the annual tropical cyclone count time series. In the view of the literature survey presented above, we consider it to be of significant importance to have an insight into the temporal behavior of the frequency of tropical cyclones. It is apprehended that outcomes of such study may help in understanding the presence of any kind of serial dependence within the time series and may generate some inputs for future predictive models. The rest of the paper is organized as follows: In Sect. 2, we have presented the data and methodology. The fitting of autoregressive models has been demonstrated in Sect. 3. In Sect. 4, we have presented a comparative study between autoregression and autoregressive neural network models. We have concluded in Sect. 5.

Data and methodology
In this section, we are going to develop the Markov chain models for the frequency of tropical cyclones over the Arabian Sea Basin, Bay of Bengal Basin, and land. The data has been collected from Cyclone eAtlas by the India Meteorological Department for the period 1891-2019, i.e., 129 years. We have categorized our study into two intensity levels: (i) CS and SCS and (ii) D, CS, and SCS, where D represents depressions having wind speed of 17 knots or more, CS represents cyclonic storms having wind speed of 34 knots or more, and SCS represents severe cyclonic storms having wind speeds of 48 knots or more. The data are converted to binary form and hence checked for two-state Markovian dependence using the χ 2 test. The order of the Markov chain implies the number of previous states influencing the current state. If X t−n , X t−(n−1) ,…., X t−2 , X t−1 influence X t , then we consider it to be a Markov chain of order n. Whenever we investigate for Markovian dependence, it is necessary to decide the value of n that gives the order of the Markov chain. Different orders may be accepted by the hypothesis testing carried out on a two-state process. However, we need to fix up the most appropriate order for further study. The order of the Markov chain (MC) is decided using the minimization of the Bayesian Information Criterion (BIC) (Ray et al. 2021;García et al. 2018). The implementation procedure is detailed in the following subsections.

Category I
In this subsection, we are going to examine the frequency of tropical cyclones considered under category I; i.e., both CS and SCS are taken into consideration. Total 668 cases are observed under this category for the study period under consideration. Each data point (x) corresponds to the total frequency of CS + SCS in 1 year within the period of study. When averaged over the entire study period, the mean frequency (x) comes out to be 5.178. Now, we convert the data series to a binary series using the following definition of a random variable: X t = 1, if x t is greater than or equal to the mean frequency and = 0 otherwise. Now, we apply the Markovian approach to this dichotomous time series to test for Markovian dependence. To do the same, we take the null hypothesis H 0 : The data are serially independent against the alternative hypothesis H 1 : There is a first-order serial dependence. The following contingency table representing the number of 4 types of transition within the dichotomous time series is generated with N ij representing the transition count from state i to state j where both i and j realize two values, i.e., 0 and 1.
Using Table 1, the transition probabilities are computed in Table 2.
In Table 2, p ij represents the transitional probability from state i to state j. The stationary probability is given by: From Table 2, using Eq. 1, we can compute the stationary probability as π 1 = 0.457. Since p 01 < π 1 < p 11 , positive serial correlation exists. Furthermore, we compute the persistence parameter r 1 = p 11− p 01 = 0.322 ≠ 0. The above computation shows that the time series is expected to have a serial correlation. To further establish the above fact and to check for first-order Markovian dependence, we carry out the χ 2 test based on the null hypothesis presented above. In this particular case, the χ 2 is computed using the formula: where n ij is the transition count as already explained and e ij = i th rowtotal×j th columntotal Totalf requency For the binary time series under consideration, χ 2 = 13.206 with degrees of freedom v = 1 and the null hypothesis is not accepted at 5% level.
Now we go to test for the second-order Markovian dependence. To do the same, the transition counts and transition probabilities are computed and presented in the Tables 3 and 4 respectively.
Since it has been observed that all the four orders of Markovian dependence are acceptable, we need to decide the best representative order of the Markov chain. For this purpose, we calculate BIC for all the orders.
where L m is the log-likelihood for order m.
From Table 6, it is apparent that BIC is minimized for the first-order Markovian dependence and hence, the first-order two-state (FOTS) model of MC is the best representative for the CS + SCS time series converted to a binary time series.

Category II
In this subsection, we apply a similar procedure as in category I. In this case, each data point represents the total frequency of D + CS + SCS in a given year. In the present case, the mean of all frequencies is 12.295. In this case, the available time series is converted to a binary time series using the following definition of a random variable: Y t = 1, if Y t is greater than or equal to the mean frequency and = 0 otherwise.
Hence, the stationary probability using Eq. 1 is computed as π 1 = 0.453. Since p 01 < π 1 < p 11 , positive serial correlation exists. We also compute the persistence parameter r 1 = p 11 − p 01 = 0.513 ≠ 0. The above computation also shows that the time series is expected to have a serial correlation. Furthermore, we calculate the Markovian dependence up to the fourth order; the χ 2 values are computed and presented in Table 7.
Since it has been observed that all the four orders of Markovian dependence are acceptable, we need to decide Table 3 Contingency table for observed transition count for secondorder Markovian dependence X t+1 = 0 X t+1 = 1 Total X t = 00 n 000 = 36 n 001 = 13 n 00. = 49 X t = 01 n 010 = 9 n 011 = 13 n 01. = 22 X t = 10 n 100 = 13 n 101 = 8 n 10. = 21 X t = 11 n 110 = 12 n 111 = 23 n 11. = 35 n .0 = 70 n .1 = 57   the best representative order of the Markov chain. Hence, BIC for all the orders is calculated again using Eq. 3. From Table 8, it is apparent that BIC is minimized for the second-order Markovian dependence and hence, the second-order two-state (SOTS) model of MC is the best representative for the D + CS + SCS time series converted to a binary time series.
The column graph plotting the BIC values for the corresponding order of the Markov chain for Category I and Category II is shown in Fig. 1.

Fitting autoregressive models
In the previous section, we have demonstrated Markov chain models for C + SCS and D + CS + SCS. In the case of CS + SCS, we have observed that the time series is characterized by a first-order Markov chain model. Also, for D + CS + SCS we have observed the second-order Markov chain to characterize the time series. Now we divide the dataset having 129 data points into a training set having 96 data points from the years 1891 to 1986 and a test set having 33 data points from the years 1987 to 2019. Based on the outcomes presented in the previous section, we apply an autoregressive approach to the time series under consideration; the time series characterized by the first-order Markov chain leads us to interpret that the state at a given time point depends on the immediately previous time point and not on the long way it has traversed to reach up to that state. Similarly, the second-order Markovian process represents a scenario where the state at a given time point depends upon the two immediately previous time points. The general autoregressive process of order K can be mathematically presented as: where µ is the mean of the time series, ϕ is the autoregressive parameter, and ϵ t+1 is a random shock or innovation which has µ ϵ = 0 and variance σ ϵ 2 and corresponds to the residual in ordinary regression. The predictand x t+1 is the value of the time series at time t + 1, and the predictor is the current value of the time series x t (Wilks, 1995). Using the training dataset, we calculate the autoregressive coefficients and with the help of the test dataset we check the goodness of fit of the abovementioned models using Willmott's index (Willmott et al. 2012)  (4)  From Table 9, it can be clearly stated that there is strong agreeability of the AR(5) model with the observed value of frequency of CS + SCS. For this AR(5) model, the mean of the residuals µ ϵ is found to be 7.81E − 16 which is approximately equal to 0. Using the Ljung-Box test, the randomness of the residuals is tested. The null hypothesis is taken that the residual is independently distributed. The test statistic is calculated using the formula: where n is the sample size, r j is the sample autocorrelation at lag j, and m is the number of lags being tested. The test statistic for sample size 96 for lag 20 is calculated to be 8.384 with a corresponding p value of 0.936 > 0.05. Thus, we do not reject the null hypothesis and conclude that the residual is independently and identically distributed; i.e., white noise  is present. Thus, the white noise variance σ ϵ 2 is calculated to be 3.034.
The column graph plotting Willmott's index of orders 1 and 2 for the above-stated autoregressive models for CS + SCS is shown in Fig. 2.
The observed and predicted values of the AR(5) model for CS + SCS are plotted in Fig. 3.
From Table 10, we observe that Willmott's index of the AR(p) models is low. Thus, these models have a lack of fit and are not good for prediction. Hence, we use a neural network approach to fit an AR(p) model to the dataset of D + CS + SCS.

Comparison of the AR model with a non-linear predictive methodology
An artificial neural network was designed to mimic a human brain. An ANN consists of an input layer consisting of a set of neurons, several hidden layers, and an output layer where each neuron is activated using an activation function. Furthermore, the output is backpropagated to the input to minimize the error function (Haykin 2009). The structure of an ANN with two hidden layers is shown in Fig. 4. In the late 1990s, ANN is known to be first implemented in the field of atmospheric sciences (Gardner and Doring, 1998). Prediction of precipitation using an ANN was first demonstrated by Hall Jun et al. (1999). Chattopadhyay (2007) has shown an ANN-based prediction of the Indian summer monsoon. Later, Azadi and Sepaskhah (2012) have also predicted the annual precipitation in the provinces of Iran using ANN. ANN is also used to model the tropospheric ozone concentration Dorling 1996, Chakraborty andChattopadhyay, 2021). Studies have been conducted to predict the air temperature with historical data considering the time lags using an ANN (Chaudhuri et al. 2005, Chattopadhyay et al. 2011). Chaudhuri et al. (2015 also showed that a multilinear perceptron can be used to nowcast visibility during fog 6 h ahead of time. Later, Nath et al. (2016) have utilized three different ANN techniques to predict the tropical cyclone activity during the post-monsoon months in the NIO. Kotal et al. (2008) have also developed a new scheme using ANN to predict the cyclone intensity for the BOB. Sharma et al. (2013) have developed the same for the Western Pacific Ocean.
In this section, a non-linear predictive model is designed using the univariate time series of frequency of cyclone for D + CS + SCS. An Autoregressive Neural Network (AR-NN) is designed for different lags. We divide the data into 75% training dataset and 25% test dataset. The input matrix of an AR-NN of order p consists of p columns for p previous states and one column for the current state. The input matrix is fed into the AR-NN and activated using the logistic function to give the predicted results for the test dataset. The logistic activation function is given by: The AR-NN models are designed up to lag 5 and Willmott's index is computed for each case. Furthermore, we compute the prediction yield (see Fig. 6) of both the AR model and AR-NN model for 5% error, 10% error, and 15% error for each lag. The prediction yield is given by the formula:  Table 11.
Totalno.ofcaseswithinx% error Totalno.of testcases The column graph comparing the AR models and AR-NN models with their Willmott's index for D + CS + SCS is shown in Fig. 5.
The column graph comparing the prediction yield of different AR models and AR-NN models for D + CS + SCS is shown in Fig. 6.

Conclusion
In the rigorous study presented in the previous sections, we have reported a Markov chain model and univariate prediction of tropical cyclones over the Arabian Sea Basin, Bay of Bengal Basin, and land collected for the period 1891-2019. The data have been categorized as per their intensity levels into two categories: (i) CS + SCS and (ii) D + CS + SCS. Here D stands for depressions with wind speed greater than or equal to 17 knots. Further details have been presented in Sect. 2 of the present paper. Since the data corresponds to continuous random variables, we have discretized them for the application of the Markov chain, which is a discrete time series approach. While discretizing, we have obtained a dichotomous time series. This methodology has been adapted to each category mentioned above and accordingly, the transition probabilities have been computed for each category to find the stationary probability. This computation has been presented in Tables 2 and Table 4. It has been observed that in each case of CS + SCS, p 01 < π 1 < p 11 and hence, it has been interpreted that a positive serial correlation exists. This has further been supplemented by a non-zero persistence parameter. To further consolidate the outcomes, we have carried out the χ 2 test to check for first-order two-state Markovian dependance. We have also carried out a similar test for Markovian dependence for the second, third, and fourth order respectively and in each case, it has been observed that the computed χ 2 is exceeding the corresponding tabular value with appropriate degrees of freedom and 5% level of significance. The results are displayed in Table 5 where it has been clearly shown that for every competing order of the two-state Markov chain model, the null hypothesis H 0 assuming serial independence is not accepted at 5% level. To choose the best order of Markov chain among the 4 orders, we have implemented the BIC minimization procedure and the results for CS + SCS are presented in Table 6. This table shows that BIC gets its minimum for the FOTS model of the Markov chain and hence, FOTS is considered to be the best representative Markovian process for CS + SCS. The outcomes are presented in Table 6. In a similar computational procedure, when carried out for D + CS + SCS, Table 7 shows that the null hypothesis of serial independence is not acceptable at the 5% level (see Table 7). However, in the case of D + CS + SCS, the BIC minimization establishes SOTS as the best representative Markov chain model (see Table 8).
To have a comparative view, we have depicted the results in Fig. 1. In the subsequent phase of the study, we have developed autoregressive models for the two categories under consideration. In Table 9, we have presented the autoregressive coefficients for CS + SCS and it is clear from this table that we have checked for autoregressive processes up to order 5. To check for the goodness of fit of the autoregressive model of a given order, we have computed Willmott's indices of orders 1 and 2. The values of Willmott's indices are also presented in Table 9 and it is visible that for AR(5) both Willmott's indices are above 0.5 and the second-order Willmott's index is above 0.8. This strongly leads us to conclude that AR (5) is the best fit autoregressive model for CS + SCS. To test the randomness of the residuals, we have implemented the Ljung-Box test where it is has been observed that the residuals are independently and identically distributed and as a consequence, we have concluded the presence of white noise within this process. We have also pictorially represented the values of Willmott's index for CS + SCS in Fig. 2. Also in Fig. 3, we have displayed the observed and predicted frequencies in the test cases, which shows that there is a significant degree of closeness in the patterns of CS + SCS frequencies in the observed and predicted cases. Contrary to what happened in CS + SCS, the autoregressive models could not perform so efficiently in the case of D + CS + SCS. Table 10 shows that an increase in the order of autoregression above 2 has resulted in decay in the values of Willmott's index. Although AR(1) and AR(2) have second-order Willmott's index close to 0.5, they cannot be interpreted as a good univariate predictive model. Considering the failure of conventional autoregressive procedure in the case of D + CS + SCS, we have developed autoregressive neural networks (AR-NN) for D + CS + SCS. It has been observed that applying logistic activation function for AR-NN we have observed that for the test cases Willmott's index is not getting any significant improvement from the conventional autoregression process. However, we have observed some significant hikes in prediction yield for the first and second orders for an acceptable prediction error of 10%. In that sense, for D + CS + SCS, implementation of the neuro-computing methodology in an autoregressive manner has given some advantage over the conventional autoregression procedure. Considering the prediction yield associated with 5% error, we have observed a significant hike in the case of the fourth-order AR-NN. Hence, in general we can say that AR-NN is a better predictive tool than conventional AR in the case of D + CS + SCS. While concluding, we would like to note that neither CS + SCS nor D + CS + SCS are characterized by serial independence. This means that somehow the frequency of cyclonic storms, severe cyclonic storms, and depressions of a given year has some influence on the subsequent year. As we incorporate depressions, the prediction of frequencies becomes more difficult. This indicates the incorporation of some degree of complexity to the system by the depressions that did not develop into cyclonic storms or severe cyclonic storms. Given the same, we propose to carry a study on the fractal behavior of the time series as our future study.