The Gegenbauer autoregressive moving-average (GARMA) model is pivotal for addressing non-additivity, non-normality, and heteroscedasticity in real-world time-series data. While primarily recognized for its efficacy in various domains, including the health sector for forecasting COVID-19 cases, this study aims to assess its performance using yearly sunspot data. We evaluate the GARMA model’s goodness of fit and parameter estimation specifically within the domain of sunspots. To achieve this, we introduce the random coefficient generalized autoregressive moving-average (RCGARMA) model and develop methodologies utilizing conditional least squares (CLS) and conditional weighted least squares (CWLS) estimators. Employing the ratio of mean squared errors (RMSE) criterion, we compare the efficiency of these methods using simulation data. Notably, our findings highlight the superiority of the conditional weighted least squares method over the conditional least squares method. Finally, we provide an illustrative application using two real data examples, emphasizing the significance of the GARMA model in sunspot research.