The covariance function approach with an iterative two-stage algorithm of Liu et al. (2000) was applied to estimate parameters for the Polish Black-and-White dairy population based on a sample of 338 808 test day records for milk, fat, and protein yields. A multiple trait sire model was used to estimate covariances of lactation stages. A third-order Legendre polynomial was subsequently fitted to the estimated (co)variances to derive (co)variances of random regression coefficients for both additive genetic and permanent environment effects. Daily and 305-day heritability estimates obtained are consistent with several studies which used both fixed and random regression test day models. Genetic correlations between any two days in milk (DIM) of the same lactation as well as genetic correlations between the same DIM of two lactations were within a biologically acceptable range. It was shown that the applied estimation procedure can utilise very large data sets and give plausible estimates of (co)variance components.