The relationship between changes of seismicity rate and static Coulomb stress is investigated at different sub–regions of the broader Aegean area. Seismic activity is studied at specific areas of interest, which are characterized by intense seismic activity and strong earthquake occurrence, known from both historical and instrumental data. The division of the area is based upon seismotectonic criteria, considering the regional kinematic properties, local seismicity and the available fault plane solutions. Coseismic stress changes are modeled and along with the tectonic loading are taken into consideration for stress change calculations. Data used for modeling contain events with magnitude greater or equal to a threshold magnitude, Mc, separately identified for each sub–region and time period. Simulations are done considering either the influence of aftershocks or declustered data. The spatial distribution of seismicity is translated into earthquake probability for both the observed and expected seismicity rates, by the application of a probability density function (PDF). Statistical process requires a normal grid superimposed on the study area. Spatial variable model parameters are calculated and then are linearly interpolated at the center of each cell of the normal grid. The dimensions of the cells are chosen in regard with the epicentral location error and the size of the catalog, such that a sufficient number of events being present in each cell and a realistic estimate of seismicity rate is done. Different values of model parameters are tried, with their limits being defined by physical and observational constraints, in order to test the sensitivity of the model in their fluctuation. Values for which results are closer to the observations are considered to express better the physical conditions and processes of the regional tectonic regime. Qualitative and quantitative correlation between the observed and the expected seismicity rates provide a test for the validity and sufficiency of the model. Qualitative correlation is done by comparing the mapped patterns of observed and expected seismicity rates, while quantitative correlation is calculated by the application of statistical tests, i.e. calculation of the correlation coefficient, r. Once the model is tested in previous cases, an estimation of the expected number of small earthquakes or the probability of a large shock to occur in the future is performed for each one of the studied sub-regions. An earthquake forecast for shocks with magnitude greater than or equal to a minimum magnitude M, is attempted for specific regions and for a settled time period, contributing to a more reliable time–dependent seismic hazard assessment.