function [pSurv,cumSurvProb] = MCSurvivorSimulation(s0,alpha,q,r,sigma,nOfYears,knockOutB,npathsIn)
% function to generate asset paths using MC simulations.

% Intial parameters

% 365 is one year: this number is hard coded.
% s0=100; % initial price
% alpha = 0.2;
% q=1.4;
% r=0.01; % IR per year
% sigma=0.3; % volatility per year
% nOfYears = 3;  % number of years
% knockOutB = 5; % price level of knock-out barrier

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if(q <= 1)
   error('!!!! ENTER: q>1 !!!!')
end

texp=1*nOfYears; % total number of years

% initialization

steps = 365*nOfYears;       % total number of days, if smallest interval time dt = 1 day.
dt = texp/steps;            % smallest time interval.

npaths = npathsIn;          % number of paths: 1000 normally used
npathg=npaths;              % store this variable
s=zeros(steps,npaths); 
s_now=zeros(npaths,1)+s0;

omega=zeros(npaths,1);
defcount=0;
randn('state',sum(100*clock));
t=dt:dt:texp;
i = 0;
for run=1:steps,
    dz=randn(npaths,1)*(dt)^(0.5);
    pomega=p(omega,t(run),q);
    
    sigmat = sigma*(t(run)^((1-q)/(2*(3-q))));
      
    s_now=s_now+r.*s_now.*dt+sigmat*s0^(1-alpha).*s_now.^(alpha).*pomega.^((1-q)/2).*dz; %TSALLIS-q-alpha PRICE EVOLUTION
    
         
    s(run,:)=s_now';
    omega=omega+pomega.^((1-q)/2).*dz;
    [n]=find(s(run,:)<knockOutB);                        %donīt continue defaulted paths!
    if ~isempty(n),  % not empty
    
        s(:,n)=[];
       
        
        omega(n)=[];
        npaths=npaths-length(n);
        defcount=defcount+length(n);
    end;
    
    if(mod(run,365) == 0)
        i=i+1;
        nDefInYear(i) = defcount;
    end
    
    s_now=s(run,:)';
end;

% Calculating the survivor probabilities per year

% defcount  % DEBUG variable: product of pSurv = 1-(defcout/npathg)

defPerYear = [nDefInYear(1) diff(nDefInYear)];

aux = npathg-nDefInYear;

totalPaths = [npathg aux(1:end-1)];

pD = defPerYear./totalPaths;

pSurv = 1-pD;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

cumSurvProb = prod(pSurv);