The integration of site effects into Probabilistic Seismic Hazard Assessment (PSHA) is a constant subject of discussion within the seismic hazard community due to its high impact on hazard estimates. To include this effect in PSHA, different approaches have been proposed by different authors. The aim of this research is to compare the hazard curves (HC) obtained for a soft, nonlinear site with two different fully probabilistic site specific seismic hazard methods: 1) The Analitical Approximation of the Full Convolution Method (AM) (Bazzurro and Cornell 2004a,b) and 2) what we call the Full Probabilistic Stochastic Method (SM). The AM computes the site-specific HC by convolving the site-specific bedrock HC, Sar(f), with a simplified representation of the probability distribution of the amplification function, AF(f) at the considered site. While the SM is nothing else that the HC 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). For this purpose, we generate sythetic hazard-consistent earthquake catalogue, apply host-to-target corrections and calculate synthetic time histories. We propagate the rock stochastic time histories, from depth to surface using two different 1D site response analysis, a linear equivalent and nonlinear codes, to evaluate the code-to-code variability. Lastly, we compute the site-specific HC with both approaches. The results found with both methods are in relatively satisfactory agreement whatever the site response code along all the studied periods. Finally, we encourage the use of stochastic simulations to integrate site effects into PSHA, since it allows to better investigate the variability of the site response better linked to more physical models, something that currently is not possible with real data due to its scarcity (specially at high acceleration levels).