Several approaches have been proposed to integrate site effects in Probabilistic Seismic Hazard Assessment (PSHA), varying from deterministic, to hybrid (probabilistic-deterministic), and finally fully probabilistic approaches. The present study compares hazard curves obtained for a soft, non-linear site with two different, fully probabilistic site specific seismic hazard methods: 1) The Full Convolution Analytical Method (AM) (Bazzurro and Cornell 2004a,b) and 2) what we call the Full Probabilistic Stochastic Method (SM). The AM computes the site-specific hazard by convolving the site-specific bedrock hazard curve, Sar(f), with a simplified representation of the probability distribution of the amplification function, AF(f) at the considered site, while the SM is built from stochastic time histories on soil corresponding to a representative, long enough catalogue of seismic events. This comparison is performed for the example case of the Euroseistest site near Thessaloniki (Greece). We generate a hazard-consistent synthetic earthquake catalogue, apply host-to-target corrections, calculate synthetic time histories with the stochastic point source approach, and scale them using an adhoc frequency dependent correction factor to fit the specific rock target hazard. We then propagate the rock stochastic time histories, from depth to surface using two different 1D site response analysis, a linear equivalent (LE) and non-linear (NL) codes, to evaluate the code-to-code variability. Lastly, we compute the probability distribution of the non-linear site amplification function, AF(f), for both site response approaches, and derive the site-specific hazard curve with both AM and SM approaches. Results are found in relatively satisfactory agreement whatever the site response code along all the studied periods. The code-to-code variability (EL and NL) is found significant, providing a much larger contribution to the hazard estimate uncertainty, than the method-to-method variability (AM and SM). However, the AM approach presents a numerical limitation, that is not encountered with the SM approach, though with a much higher computational price. The use of stochastic simulations to integrate site effects into PSHA allows to better investigate the variability of the site response physics, and a good parameterization of the input parameters, something that currently is not possible with real data due to its scarcity especially at high acceleration levels.