In many real world situations there is no reason to believe that the time series observations are normally distributed. Therefore the estimation of the error distribution, estimation of the parameters of the process and prediction for non-Gaussian autoregressive time series models are of importance. In our development, it is assumed that there exists a transformation that transforms the error distribution to the normal distribution, and that this transformation can be represented by a regression spline function. The transformation, and hence, the error distribution is estimated by an estimation procedure based on the spline regression of the residual quantiles on the corresponding normal scores. The limiting distribution of the estimator of the vector of spline parameters is derived. Given an estimator of the error distribution, a nonlinear maximum likelihood estimation procedure is used to obtain an improved estimator of the vector of autoregressive parameters and to construct confidence intervals for predictions. The limiting distribution of the estimator of the vector of autoregressive parameters is shown to be normal;Monte Carlo simulations conducted with different error distributions demonstrate that the procedure performs well for finite samples. The estimated quantiles of the error distributions are close to the true quantiles, although there are significant biases. The Monte Carlo study indicates that the variance of the autoregressive parameter estimator for the spline method is smaller than the variance of the corresponding least squares estimator for skewed and long tailed error distributions such as chi-square distributions. For normal errors the spline method is only moderately inferior to the least squares method. In the Monte Carlo study confidence intervals for predictions were constructed using normal distribution theory and using the distribution estimated by the spline method. For chi-square and mixture of normals error distributions, coverages of the intervals based on the spline method are superior to the coverages of the intervals based on the normal distribution. For normal errors the performance of the spline intervals is very close to that of intervals based on the normal distribution.