Estimating the arrival rate function of a non-homogeneous Poisson process based on observed arrival data is a problem naturally arising in many applications. Cubic spline functions are particularly well-suited to represent such estimates since they are sufficiently versatile and easy to handle computationally. In this paper we present an optimization model to obtain cubic spline estimations based on the maximum likelihood principle. An important feature of any such model is that in order for us to be able to interpret the obtained splines as arrival rate functions they have to be nonnegative over the time interval of interest. We ensure this using a result of Nesterov characterizing nonnegative polynomials via positive semidefinite matrices. We also describe versions of our model allowing for periodic arrival rate functions and input data of limited precision. Numerical results based both on real-life and artificially generated data sets are also presented. 1