Abstract Temporal and spatial variation in phenotypic selection due to changing environmental conditions is of great interest to evolutionary biologists, but few existing methods estimating its magnitude take into account the temporal autocorrelation. We use state-space models (SSMs) to analyse phenotypic selection processes that cannot be observed directly and use Template Model builder (TMB), an R package for computing and maximizing the Laplace approximation of the marginal likelihood for SSM and other complex, nonlinear latent variable model via automatic differentiation. Using a long-term great tit (Parus major) dataset, we fit several SSMs and conduct model selection based on Akaike information criterion (AIC) to assess the support for fluctuated directional or autocorrelated stabilizing selection on breeding time of the great tit population. Our results show that there is directional selection on the probability of breeding failure, and stabilizing selection on the mean number of fledglings. This selection for early laying date is consistent with a previous study of the same population. We also estimate the variation and autocorrelation in other parameters of the fitness functions, including the width and height, and found the height and location of annual fitness function are autocorrelated with significant variation, while the width can be assumed to constant over time. Using TMB to fit SSMs, we are able to estimate additional parameters compared to previous methods, all without requiring a substantial increase in computational resources. Furthermore, our specification of complex nonlinear model structure benefits greatly from the flexibility of model formulation with TMB. Therefore, our approach could be directly applied to estimating even more complicated phenotypic selection processes induced by environmental change for other species.