P>We present an algorithm for automatic P-phase arrival time determination for local and regional seismic events based on higher order statistics (HOS). Using skewness or kurtosis a characteristic function is determined to which a new iterative picking algorithm is applied. For P-phase identification we apply the Akaike Information Criterion to the characteristic function, while for a precise determination of the P-phase arrival time a pragmatic picking algorithm is applied to a recalculated characteristic function. In addition, an automatic quality estimate is obtained, based on the slope and the signal-to-noise ratio, both calculated from the characteristic function. To get rid of erroneous picks, a Jackknife procedure and an envelope function analysis is used. The algorithm is applied to a large data set with very heterogeneous qualities of P-onsets acquired by a temporary, regional seismic network of the EGELADOS-project in the southern Aegean. The reliability and robustness of the proposed algorithm is tested by comparing more than 3000 manually derived P readings, serving as reference picks, with the corresponding automatically estimated P-wave arrival times. We find an average deviation from the reference picks of 0.26 +/- 0.64 s when using kurtosis and 0.38 +/- 0.75 s when using skewness. If automatically as excellent classified picks are considered only, the average difference from the reference picks is 0.07 +/- 0.31 s and 0.07 +/- 0.41 s, respectively. However, substantially more P-arrival times are determined when using kurtosis, indicating that the characteristic function derived from kurtosis estimation is to be preferred. Since the characteristic function is calculated recursively, the algorithm is very fast and hence suited for earthquake early warning purposes. Furthermore, a comparative study with automatically derived P-readings using Allen's and Baer & Kradolfer's picking algorithms applied to the same data set demonstrates better quantitative and qualitative performance of the HOS approach. This study shows, that precise automatic P-onset determination is feasible, even when using data sets with very heterogeneous signal-to-noise ratio.