Approach of the value of a rent when non-central moments of the capitalization factor are known: an R application with interest rates following normal and beta distributions Ratio Mathematica Volume 41, 2021, pp. 119-136 119 A unified shape model for sunspot number cycles Beena Girija Puthumana* Sabarinath Amaranathan† Anilkumar Ajimandiram Krishnankuttynair‡ Abstract We proposed a model which can unify many of the shape models existing in the literature and show that the shape of the sunspot number cycle can be described as a product of a polynomial and a negative exponential function. The proposed model has certain free parameters, which need to be estimated from the observed sunspot number data. Since all the models reviewed in this paper are a product of a polynomial and a negative exponential along with a number of parameters, we have seen that all these models can be derived from a modified generalized Gamma probability density function by transforming certain parameters and fixing certain parameters. In this paper, we estimate the parameters of the model from the revised monthly averaged sunspot numbers available in the SIDC website. Finally, a preliminary level prediction has also been attempted to forecast the characteristics of sunspot number cycle 25. Keywords: sunspot numbers; gamma distribution; free parameter; etc. 2010 AMS subject classification: 70F15; 97M10.§ *Assistant Professor, St. Gregorios College, Kottarakkara, Kerala, India; beenamabhi@gmail.com. †Scientist, Vikram Sarabhai Space Centre, Thiruvananthapuram, Kerala, India; a_sabarinath@yahoo.co.in. ‡ Director, DSSAM, ISRO Head Quarters, Bangalore, Karnataka, India; ak_anilkumar@isro.gov.in. § Received on xxx, xxx. Accepted on xxx, xxx. Published on xxx, xxx. doi: xxx. ISSN: xxx. eISSN: doi. 10.23755/rm.v41i0.671.. ©The Authors. This paper is published under the CC-BY licence agreement. Beena G P, Sabarinath A, Anilkumar A K 120 1. Introduction Solar activity is the key factor which drives the space weather. Refined modelling and accurate prediction of the solar activity intensity has been an important activity of space faring agencies across the world due to the impact of solar activity on satellites as well as on weather (Haigh, 2007, Hathaway et.al., 2004). Solar flux causes the upper atmosphere density variation and in turn it affects directly the lifetime of the Earth-orbiting satellites especially in the low Earth orbit. Solar activity intensity has been measured as the number of dark spots, called sunspot numbers, appears in the visible solar disc through direct observation since 1749 onwards. Irrespective of the measurement interval (daily, monthly, and yearly), a definite pattern is existing in the sunspot number time series. Accurate predictions of the intensity of solar activity are increasingly important as we become more reliant upon satellites in low-Earth orbits, which provide crucial contribution in communication, national defence and Earth mapping. Also, such satellites provide an abundance of scientific data. During higher solar activity period, the increased ultraviolet emission from Sun heats up the Earth’s upper atmosphere and this causes the atmosphere to expand and results in the increased drag on low Earth orbits satellites, thereby leading to early decay into the Earth’s atmosphere. Therefore, better predictions of solar activity are essential to help mission planning and design of satellites (Vallado et.al., 2014). Sunspot number cycle time series is one of the longest time series which was studied by many experts. First of all, this time series is non stationary, cyclic and highly nonlinear in the time domain. The more interesting and difficult part to deal with is the high dispersion between consecutive observations (Withbroe,1989) which in fact makes the prediction of sunspot numbers tedious. Many attempts to model and predict the future behaviour of the solar activity are well documented in the literature. Depending on the nature of the prediction methods, we can classify the methodology in to five classes as: 1) Curve fitting 2) Precursor 3) Spectral 4) Neural Networks and 5) Climatology. The first attempt using the curve fitting methodology was by the McNish- Lincolon curve fitting (Hathaway, D. H, 2015). Subsequently, several authors have studied the highly nonlinear behaviour of sunspot number cycle and proposed various models to handle the studies related to the prediction. Many mathematical functions have been appeared in the literature as model of the shape of the sunspot number cycle. Due to the exponential rise and decay of sunspot number cycle, a model involving exponential function was proposed by Nordemann (1992). The bell shape and the asymmetry along the peak amplitude of most of the sunspot number cycle were explored and there by a suitable mathematical function was introduced by Hathaway et.al. (1994). Few A Unified Shape Model for Sunspot cycles 121 statistical probability distribution functions were also proposed by various authors (Sabarinath et.al.,2008, Du et.al., 2011, Li et.al., 2017, Sabarinath et.al., 2020) for modeling the shape of sunspot number cycles. De Meyer (1981) proposed a model using periodic functions. As far as prediction of a future sunspot number cycle is concerned, statistical averaged models are used as an initial estimate of the future sunspot number cycle. 2. Existing Models Several authors developed different mathematical functions to describe the shape of the sunspot number cycle. In particular Stewart and Panofsky (SP) (1938) proposed a function for the shape of a sunspot number cycle with the form, 𝑅(𝑡) = 𝑎(𝑡 − 𝑡0) 𝑏 𝑒−𝑐(𝑡−𝑡0) (1) where 𝑎, 𝑏, 𝑐 and 𝑡0 are parameters that vary from cycle to cycle and 𝑡 is the independent variable represents time. The important thing to be noticed is that, this model resembles a power law for the rising phase of a sunspot number cycle and an exponential for the declining phases of a cycle. Nordemann (N) (1992) proposed another fit. He used the solution of the differential equation 𝑑𝑁 𝑑𝑡 = 𝐾𝑁, in analogy with the nuclear decay process. The declining phase (maximum to minimum) of a sunspot number cycle is represented by 𝑁 = 𝑁0𝑒 𝐾𝑡 , 𝐾 < 0 (2) and the solution of 𝑑𝑁 𝑑𝑡 = 𝐴 + 𝐾𝑁, is used to represent the first phase or ascent phase (minimum to maximum) of a sunspot number cycle. Thus, the model for the ascent phase is: 𝑁 = 𝐴 𝐾 (1 − 𝑒𝐾𝑡 ), 𝐾 < 0 (3) where 𝑁 represents sunspot numbers, 𝐾 is the decay constant and 𝐴 is a production parameter. Hathaway et.al. (H) (1994) established a model with four parameters along with a measure for the goodness of fit. The functional form is: 𝑓(𝑡) = 𝑎(𝑡 − 𝑡0) 3 𝑒 [ (𝑡−𝑡0) 2 𝑏2 ] − 𝑐 (4) where, 𝑎 represents the amplitude; 𝑏 represents the time in months and 𝑡0 denotes the starting time; 𝑐 gives the asymmetry of the sunspot number cycle. This function is derived from Stewart and Panofsky model, but requires two Beena G P, Sabarinath A, Anilkumar A K 122 parameters for each sunspot number cycle. Where 𝑐 = 0.71, and 𝑏 is a dependent parameter given by, 𝑏 = 27.12 + 25.15 [𝑎 × 103] 1 4 (5) Sabarinath et.al. (S) (2008) used a binary mixture of a modified Laplace Distribution. Laplace Distribution is a function 𝑓 of two parameters 𝑚 and 𝑠 is given by, 𝑓(𝑥) = 1 2𝑠 𝑒 −|𝑥−𝑚| 𝑠 (6) where, 𝑚 is the location parameter and 𝑠 is the scale parameter. They modified this form and used the binary mixture of this distribution and reducing the number of parameters into six (later reduced to two floating parameters) to fit the predominant double peaks during the high solar activity regime of a sunspot number cycle. Modified final model is: Volobuev (V) (2009) found a function similar to that used by Stewart and Panofsky to fit the shape of sunspot number cycle that requires only one parameter for each cycle. The empirical model used is: 𝑅 = ( 𝑡 − 𝑡0 𝑇𝑠 ) 2 𝑒 −( 𝑡−𝑡0 𝑇𝑑 ) 2 (8) We can see that this model is also similar to that of Stewart and Panofsky (1938) by putting 𝑏 = 2 and modifying the growth multiplier and decay multiplier properly by introducing the new parameters 𝑇𝑠 and 𝑇𝑑. Du (2011) suggested modified Gaussian function with four parameters viz. peak size 𝐴, peak timing 𝑡𝑚, width 𝐵, and asymmetry 𝛼, in the form: 𝑅(𝑡) = 𝐴 𝑒𝑥𝑝 ( −(𝑡 − 𝑡𝑚) 2 2𝐵2[1 + 𝛼(𝑡 − 𝑡𝑚)] 2 ) (9) Li et al (L) (2017) used a binary mixture of Gaussian function (suggested by Du), 𝑓(𝑥) = 𝐴1 𝑒𝑥𝑝 (− (𝑥 − 𝑚1) 2 𝑠1 ) + 𝐴2 𝑒𝑥𝑝 (− (𝑥 − 𝑚2) 2 𝑠2 ) (10) this model has six parameters. Peak sizes are denoted by and 𝐴1, 𝐴2, 𝑠1 and 𝑠2 represents gradients and peak time is represented by 𝑚1 and 𝑚2. Sabarinath et. al. (SB) (2020) fit the full sunspot number cycle perfectly with modified Maxwell Boltzmann probability distribution function. The final modified model is: 𝑓(𝑡) = 𝐴1 33.2 𝑒𝑥𝑝 ( −|𝑡 − 𝑡0 − 41.7| 16.6 ) + 𝐴2 46 𝑒𝑥𝑝 ( −|𝑡 − 𝑡0 − 67.3| 23 ) (7) A Unified Shape Model for Sunspot cycles 123 𝑓(𝑥; 𝛼; 𝐴) = 𝐴 𝛼3 √ 2 𝜋 𝑥2𝑒 − 𝑥2 2𝛼2 (11) where 𝐴 is the area parameter, 𝛼 > 0. All these models discussed so far has a common form that is these models are a product of a polynomial and a negative exponential function. Since the solar activity like process can be modelled by a bell-shaped curve viz., Gamma family of probability density function, we propose that the sunspot number cycles can be effectively modelled as a generalized Gamma distribution with certain free parameters. These free parameters can be estimated in the maximum likelihood sense from the sunspot data. Also, all these models discussed so far can be derived from this generalized Gamma distribution model as special cases. In the next section we brief on the derivation of the proposed model. 3. Methodology New model-Generalized Gamma distribution The Gamma distribution is often used to describe variables bounded on one side. A version of this distribution is obtained by adding a third parameter and gets the generalized Gamma distribution (Walck, 2001). Probability density function is, 𝑓 (𝑥; 𝑝, 𝑞, 𝑟) = 𝑝𝑟 𝛤(𝑞) (𝑝𝑥)𝑞𝑟−1𝑒−(𝑝𝑥) 𝑟 (12) where 𝑝 (a scale parameter) and 𝑞 are the real positive parameters and a third parameter 𝑟 has been added (𝑟 = 1 for the ordinary Gamma distribution) to generalize the Gamma distribution. This new parameter takes any real value but normally we consider the case where 𝑐 > 0 Put the following substitutions, 𝑐 = 𝑝𝑟 𝛤(𝑞) 𝑝𝑞𝑟−1  = 𝑞𝑟 𝐾 = 𝑝𝑟 𝛿 = 𝑟 (13) Then, Equation (12) becomes 𝑓 (𝑥; 𝑐, 𝑘, 𝛼, 𝛿) = 𝑐𝑥𝛼−1𝑒−𝑘𝑥 𝛿 (14) where 𝑘 is the scale parameter, 𝛼 is the shape parameter and 𝛿 is the location parameter. Depending on the values of the parameters we can arrive all the Beena G P, Sabarinath A, Anilkumar A K 124 existing models. Table1 list the existing models and their derivation into Equation (14). Once we fix a model, the next step is to evaluate the best estimates of the model parameters in statistical parameter estimation sense. Here our measurement data is the monthly averaged sunspot numbers. The model we indent to fit over this data is given in Equation (14). We estimate the parameters by least square method. Using simple random search technique, we estimate the parameters. The mathematical algorithms used are given in detail in the next part. Data In the present study, we use the monthly averaged sunspot numbers, and these sunspot numbers for all the 24 sunspot number cycles were used. On July 1st, 2015, the sunspot number series has been replaced by a new improved version called version 2.0 data, that includes several corrections of past in homogeneities in the sunspot number time series. In the new version 2.0 data, the conventional 0.6 Zürich scale factor, has been replaced by a factor of 1/0.6. This scale change, when combined with the recalibration, leads to a net increase of about 45% (correction variable with time) of the most recent part of the series, after 1947. This data can be obtained from https://www.bis.sidc.be/silso/DATA/SN_m_tot_V2.0.txt. Estimation techniques The function in which parameters to be estimated is, 𝑓 (𝑥; 𝑐, 𝑘, 𝛼, 𝛿) = c𝑥𝛼−1𝑒−𝑘𝑥 𝛿 (15) The maximum likelihood estimates of the parameters 𝛼, 𝛿 and 𝑘 are considered to be the best unbiased, consistent and sufficient estimate of the parameters (Sorenson, 1980). Practically, the least square estimate is considered to be the maximum likelihood estimate. The simple mathematical procedure to estimate the parameters is to minimize the sum of squared error function 𝐽. 𝐽 = ∑ 𝑒𝑟 2 𝑟 (16) where 𝑒𝑟 is the error. The minimum of 𝐽 can be found by differentiating 𝐽 with respect to the parameters 𝛼, 𝛿 and 𝑘. In the present study, if we consider without loss of generality, a sunspot number cycle having a length of 132 months (11 year), and if we assume {𝑠𝑛: 𝑛 = 1,2, … ,132} as the realised sunspot number values, then the function 𝐽 can be written as, https://www.bis.sidc.be/silso/DATA/SN_m_tot_V2.0.txt A Unified Shape Model for Sunspot cycles 125 𝐽 = ∑[𝑠𝑛 − 𝑓(𝑥𝑛, 𝛼, , 𝑘)] 2 132 𝑖=1 (17) where, 𝑥𝑛 = 1,2, … ,132 represents the months for each 𝑛 = 1,2, … . ,132. Then our objective is to compute and solve 𝛼, 𝛿 and 𝑘 from 𝜕𝐽 𝜕𝛼 = 0 (18) 𝜕𝐽 𝜕 = 0 (19) 𝜕𝐽 𝜕𝑘 = 0 (20) Analytically solving the Equations (18) to (20) for 𝛼, 𝛿 and 𝑘 is cumbersome. Hence, we go with numerical procedures for estimating the parameters. Monte Carlo based simple random search-based procedure is considered here to estimate the parameters. This procedure is described below as an algorithm. Step-1. Start with a search region 𝛼, 𝛿 and 𝑘. Let 𝑆𝛼 , 𝑆𝛿 and 𝑆𝑘 are the bounded search regions of 𝛼, 𝛿 and 𝑘 respectively. Our objective is to find an 𝛼0 ∈ 𝑆𝛼 , 𝛿0 ∈ 𝑆𝛿 and 𝑘0 ∈ 𝑆𝑘 such that, 𝐽𝛼0,𝐴0 = ∑[𝑠𝑛 − 𝑓(𝑥𝑛, 𝛼0, 𝛿0, 𝑘0)] 2 132 𝑖=1 (21) is minimum. That is, 𝐽𝛼0,𝛿0,𝑘0 ≤ 𝐽𝛼,,𝑘 (22) for any 𝛼 ∈ 𝑆𝛼, 𝛿 ∈ 𝑆𝛿 and 𝑘 ∈ 𝑆𝑘 Step-2. Start with a random initial value of 𝛼 in 𝑆𝛼, 𝛿 in 𝑆𝛿 and 𝑘 in 𝑆𝑘. Compute 𝐽 and in each iteration keep the minimum value of 𝐽, 𝛼, 𝛿 and 𝑘. After a very large number of iterations take the value of 𝛼, 𝛿 and 𝑘 corresponds to the global minimum value of 𝐽. 4. Results Estimates for the parameters Table 2 shows typical converged values of the four model parameters 𝑐, 𝛼, 𝑘 and 𝛿. If we do a Monte Carlo based estimation of these parameters, due to the initial random number variation the optimum value differ numerically due to different starting points. But the variation is insignificant. This was proved in many Monte Carlo based optimizations (Ji et.al., 2006). Hence without loss of generality we consider a typical Monte Carlo run and a converged value of the parameters for further analysis. Table 2 gives one such value. The range of Beena G P, Sabarinath A, Anilkumar A K 126 search or feasible region was found by trial-and-error method. The range considered for the simulation is given below. 0.0001 ≤ 𝑐 ≤ 0.0030, 3 ≤ 𝛼 ≤ 7, 0 ≤ 𝑘 ≤ 0.9, and 0 ≤ 𝛿 ≤ 3.1. Model Name and Year in which it is proposed Functional form Parameters in Equation (14) Generalized form 𝑐 𝛼 𝛿 𝑘 SP (1938) a(𝑡 − 𝑡0) 𝑏 𝑒𝑥𝑝(−𝑐(𝑡 − 𝑡0)) free free 1 free a𝑥 𝛼 𝑒𝑥𝑝 (−𝑐𝑥) N (1992) 𝑁 = 𝑁0𝑒𝑥𝑝(𝑘𝑡) 𝐴 𝐾 (1 − 𝑒𝑥𝑝(−𝑘𝑡)) free 1 1 free a𝑒𝑥𝑝 (−𝑐𝑥) H (1994) 𝑎(𝑡 − 𝑡0) 3 𝑒𝑥𝑝 ( (𝑡 − 𝑡0) 2 𝑏2 ) − 𝑐 free 4 2 free 𝑎𝑥3𝑒𝑥𝑝(−𝑘𝑥2) S (2008) 𝐴1 33.2 exp ( −|𝑡 − 41.7| 16.6 ) + 𝐴2 46 𝑒𝑥𝑝 ( −|𝑡−67.3| 23 ) free 1 1 fixed 𝑐1𝑒𝑥𝑝(−𝑘1x) +𝑐2𝑒𝑥𝑝(−𝑘2x) V (2009) (𝑡 − 𝑡0) 2 𝑇𝑠 2 𝑒𝑥𝑝 ( −(𝑡 − 𝑡0) 2 𝑇𝑑 2 ) free 3 2 free 𝑐𝑥2𝑒𝑥𝑝(−𝑘𝑥2) Du (2011) 𝐴 𝑒𝑥𝑝 (− (𝑡 − 𝑡𝑚) 2 2𝐵2(1 + 𝛼(𝑡 − 𝑡𝑚) 2 ) free 1 2 free 𝑐𝑒𝑥𝑝(−𝑘𝑥2) L (2017) 𝐴1𝑒𝑥𝑝 ( −(𝑡−𝑚1) 2 𝑠1 )+𝐴2𝑒𝑥𝑝 ( −(𝑡−𝑚2) 2 𝑠2 ) free 1 2 free 𝑐1𝑒𝑥𝑝(−𝑘1𝑥 2) +𝑐2𝑒𝑥𝑝(−𝑘2𝑥 2) SB (2020) 𝐴 𝛼3 √2 𝜋⁄ 𝑡 2 𝑒𝑥𝑝 ( −𝑡2 2𝛼2 ) free 3 2 free 𝑐𝑥2𝑒𝑥𝑝(−𝑘𝑥2) Table 1. Different models, their parameters and its values. Figure 1, 2, and 3 shows the model fit of the model on sunspot number cycles 13, 23 and 24, respectively. A Unified Shape Model for Sunspot cycles 127 Figure 1. Generalized Gamma distribution fit on the monthly averaged sunspot number cycle 13. Figure 2. Generalized Gamma distribution fit on the monthly averaged sunspot number cycle 23. Beena G P, Sabarinath A, Anilkumar A K 128 Figure 3. Generalized Gamma distribution fit on the monthly averaged sunspot number cycle 24. Analysis of the parameters The parameters are estimated on each of the sunspot number cycle independently. Figures 4 and 5 show these estimated model parameters for sunspot number cycle 1 to 24. The parameters are estimated in the maximum likelihood (ML) sense using the random search method. In table 2 Column number 2 to 5 gives the ML estimate of the parameters 𝑐, 𝛼, 𝑘 and 𝛿. Column 6 gives the coefficient of determination 𝑅2 value. From these values one can see that the goodness of fit of the model is fair and on all modern sunspot number cycles are of high degree, since the coefficient of determination is greater than 0.8. The important thing to be noted is that, since the model is derived from the Gamma distribution probability density function as given in Equation (12), there must be a correlation among the model parameters. As evident through the set of Equations (13) the correlation between the model parameters is given in Table 3. 𝛼 has a very high correlation with 𝑘 and 𝛿. Figure 6 and 7 shows this correlation along with the linear regression model derived out of this correlation. of course, 𝛼 and 𝑘 has a positive correlation and between 𝛼 and 𝛿 a negative correlation. The corresponding linear regression fits are given in Equation (23) and (24). A Unified Shape Model for Sunspot cycles 129 𝑘 = 0.29𝛼 − 1.3 (23) and 𝛿 = −0.32𝛼 + 2.7 (24) Now, substitute Equation (23) and (24) in Equation (14), we can reduce the proposed model into a two-parameter model. Again, if we re-estimate the two parameters in a maximum likelihood sense, we can again obtain a correlation between the parameters, there by the model reduce to a one parameter model as proposed in the shape of the sunspot number cycle-a one parameter fit by Volobuev (2009). In fact, the reduction of model parameters into a single parameter does not add any predictive power in the characterisation of a sunspot number cycle via the prediction of the peak amplitude, location of the sunspot maximum, cycle length etc. Hence, an indicative parameter for these characters of a sunspot number cycle is essential in a sunspot model. So minimum two parameters, a location parameter and a scale parameter must be there in a sunspot model. If there is a shape parameter, it can characterize the degree of asymmetry present in a cycle. Figure 4. Variation of the ML estimate of the parameter 𝑐 over different sunspot number-cycles Coefficient of determination (𝑅2) is considered as one of the measures of goodness of fit for a regression fit. For each cycle the coefficient of determination for the best fit are computed and is given in Table 2, column 6 and the pictorial version is given in Figure 8. Beena G P, Sabarinath A, Anilkumar A K 130 Sunspot number cycle Number 𝑐 α 𝑘 δ Coefficient of determination 1 0.00095 4.62004 0.03367 1.10253 0.7 2 0.00259 5.86033 0.52436 0.69427 0.7 3 0.00104 6.79747 0.80867 0.65381 0.9 4 0.00278 5.83371 0.56733 0.66284 0.9 5 0.00213 4.21489 0.01729 1.22798 0.8 6 0.00144 3.97518 0.00191 1.62772 0.6 7 0.00114 4.23167 0.00323 1.53445 0.7 8 0.00014 6.71511 0.47050 0.72610 0.8 9 0.00044 5.30111 0.10092 0.93291 0.8 10 0.00033 6.01638 0.37419 0.72687 0.9 11 0.00153 5.29585 0.14993 0.89630 0.9 12 0.00119 4.53186 0.01235 1.33943 0.8 13 0.00226 6.11178 0.74136 0.63797 0.9 14 0.00277 5.06902 0.26096 0.77246 0.7 15 0.00278 4.34677 0.00650 1.48307 0.8 16 0.00283 4.81611 0.08682 0.98842 0.8 17 0.00096 4.98223 0.04514 1.10391 0.9 18 0.00238 4.96288 0.06138 1.07036 0.9 19 0.00298 5.12209 0.08539 1.02239 0.9 20 0.00242 4.74031 0.08620 0.95019 0.9 21 0.00236 5.06018 0.08613 1.00536 0.9 22 0.00205 5.00316 0.06184 1.07074 0.9 23 0.00288 4.78601 0.07280 1.00605 0.9 24 0.00280 4.84851 0.09596 0.98271 0.8 Table 2. The estimated model parameters for sunspot number cycles 1 to 24. It may be seen that the 𝑅2 value for most of the cycles are greater than 0.8. Especially, for modern cycles (cycle 15 to 24) the 𝑅2 values are greater than 0.8. This shows that the generalized Gamma model is best for modelling the shape of a sunspot number cycle. 𝑐 𝛼 𝑘 𝛿 𝑐 1 -0.27 -0.09 -0.02 𝛼 1 0.89 -0.86 𝑘 1 -0.78 𝛿 1 Table 3. The correlation among the estimated model parameters. A Unified Shape Model for Sunspot cycles 131 Figure 5. Variation of the ML estimate of the parameter 𝛼, 𝑘, and 𝛿 over different sunspot number cycles. Figure 6. Linear correlation between the model parameters 𝛼 and 𝑘. Beena G P, Sabarinath A, Anilkumar A K 132 Figure 7. Linear correlation between the model parameters 𝛼 and 𝛿. Figure 8. Coefficient of determination (𝑅2) of the ML fit on different cycles. A Unified Shape Model for Sunspot cycles 133 Prediction of sunspot number cycle 25 As an attempt has been made to predict the shape of the sunspot number cycle 25. Before attempting the predict cycle 25 from the proposed model, a direct comparison of the model of all the sunspot number cycles from 1 to 24 has been done. Figure 9 shows this model comparison. It may be noted that, the length of the cycle is not being considered here. Without loss of generality one can assume the length as 132 months. The variation in amplitude and the location of the peak amplitude varies between cycles. The range of variation of peak amplitude is from 70 to 300 units of sunspot numbers. Since the range of variation of the peak amplitude is quite large, the problem of prediction of cycle 25 is cumbersome from the previous cycle’s model parameters alone. Hence, we make two kinds of predictions which are statistically more probable forecasts. They are (1) Since the parameters 𝛼, 𝑘 and 𝛿 has lesser variation if we consider cycle 16-24, we fix these parameters as their average over cycle 16 to 24 and re-estimated the other parameter 𝑐 alone for these cycles. Then for cycle 25, the parameter is 𝑐 taken as the average of the re-estimated values of 𝑐 for cycle 16 to 24. This value is 0.002. The remaining parameters are, 𝛼 = 4.9246, 𝑘 = 0.0757, 𝛿 = 1.02. The prediction of cycle 25 in this direction is plotted in Figure 10. Figure 9. Model of all the 24 Sunspot number cycles. This prediction shows the peak amplitude as 157 units occurring at 47 months from the beginning of sunspot number cycle 25. The second methods are (2) keeping the parameters 𝛼, 𝑘, and 𝛿 similar to cycle 24 and keeping the Beena G P, Sabarinath A, Anilkumar A K 134 parameter 𝑐 as taken in method (1). In this way the peak amplitude of cycle 25 is 81 units of sunspot numbers occurring at 44 months from the beginning of the cycle 25. These two predictions can be considered as a band, inside which the actual cycle 25 may occur. Figure 10. Prediction of Sunspot number cycle 25 Maximum peak occurring at 44 months with a sunspot number value of 81 units and 47 months with a sunspot value of 157 units. 5 Conclusion Proposed a model which can unify many of the shape models existing in the literature. Also, it is shown that the shape model of sunspot number cycle can be described by a product of a polynomial and a negative exponential function. Since all the models reviewed in this paper are a product of a polynomial and a negative exponential, we proposed that all these models can be derived from the generalized Gamma probability density function by giving suitable parameter values. In this paper, we derived the existing models from the proposed generalized Gamma model and estimated the parameters of the proposed model from the revised version-2 monthly averaged sunspot numbers available in the SIDC’s website. Prediction of sunspot number cycle 25 shows that the peak amplitude of cycle 25 can vary between 81 to 157 units of sunspot numbers and this peak amplitude may occur between 44 to 47 months from the beginning of cycle 25. In actual date, this shows cycle 25 may peak during August 2023 to November 2023. A Unified Shape Model for Sunspot cycles 135 References [1] Boslaugh, S. Watters, P.A. Statistics in a nutshell. Sebastopol, CA: OReilly Media. 2008. [2] De Meyer, F. Mathematical modelling of the sunspot cycle. Solar Physics, 70(2): 259-272, 1981. [3] Du, Z. The shape of solar cycle described by a modified Gaussian function. Solar Physics, 273(1): 231-253, 2011. [4] Hathaway, D. H. Wilson, R. M. Reichmann, E. J. The shape of the sunspot cycle. Solar Physics, 151(1): 177-190, 21, 1994. [5] Haigh, J.D. The Sun and the Earth’s climate. Living reviews in solar physics. 4(1): 1-64, 2007. [6] Hathaway, D. H. Wilson, R. M. What the sunspot record tells us about space climate. Solar Physics, 224(1): 5-19, 2004. [7] Hathaway, D. H. The solar cycle. Living reviews in solar physics, 12(1): 4, 2015. [8] Ji, M., Klinowski, J. Taboo evolutionary programming: a new method of global optimization. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 462: 3613–3627, 2006. [9] Li, F. Y. Xiang, N. B. Kong, D. F. Xie, J. L. The shape of solar cycles described by a simplified binary mixture of Gaussian functions. The Astrophysical Journal, 834(2): 192, 2017. [10] Nordemann, D. J. R. Sunspot number time series: exponential fitting and solar behavior. Solar physics, 141(1): 199-202, 1992. [11] Sabarinath, A. Anilkumar, A.K. Modeling of Sunspot numbers by a modified binary mixture of Laplace distribution functions. Solar Physics, 250(1): 183-197, 2008. [12] Sabarinath, A. Beena, G. P. Anilkumar, A. K. Modelling the shape of sunspot cycle using a modified Maxwell-Boltzmann probability distribution function. Ratio Mathematica, 39: 33-54, 2020. [13] Sorenson, H.W. Parameter Estimation-Principles and Problems. Marcel Dekker, New York. 1980. [14] Stewart, J. Q. Panofsky, H. A. A. The Mathematical characteristics of sunspot variations. The Astrophysical Journal, 88: 385, 1938. [15] Vallado, D. A. Finkleman, D. A critical assessment of satellite drag and atmospheric density modeling. Acta Astronautica, 95: 141-165, 2014. Beena G P, Sabarinath A, Anilkumar A K 136 [16] Volobuev, D. M. The shape of the sunspot cycle: A one-parameter fit. Solar physics, 258(2): 319-330, 2019. [17] Withbroe,G.L. Solar activity cycle-History and predictions. Journal of Spacecraft and Rockets, 26(6): 394-402. 22. 1989. [18] Walck, C. Hand-book on statistical distributions for experimentalists. Internal Report SUF–PFY/96–01, Particle Physics Group, University of Stockholm. 2001.