function Probability = UnconditionalProbability(Aobs,J,model,interp)
% UNCONDITIONALPROBABILITY Calculates the probability of the first observation

%% Initialization
delta       = model.params(3);
StockInterp = interp.cx{1};
nburn       = 20; % Number of burn-in periods

%% Random shocks
if isfield(model,'harvest') && ~isempty(model.harvest)
  %% Use provided matrix of random shocks
  harvest          = model.harvest;
  [nrep, nper]     = size(harvest);
else
  %% Generate a matrix of random shocks
  rng(1)
  nrep             = 1E4;
  nper             = 100+nburn;
  harvest          = NaN(nrep,nper);
  harvest(:,2:end) = model.funrand(nrep,nper-1);
end

nasymp = nrep*(nper-nburn);

%% Model simulation
S        = NaN(nrep,nper);
S(:,1)   = max(StockInterp(zeros(nrep,1)),0);
for t=2:nper
  A      = (1-delta)*S(:,t-1)+harvest(:,t);
  S(:,t) = max(StockInterp(A),0);
end

Sasymp = reshape(S(:,(nburn+1):nper),nasymp,1);

%% Probability of the first point
Probability = NaN(size(Aobs));
for i=1:length(Aobs)
  omega          = Aobs(i)-(1-delta)*Sasymp;
  ProbMC         = exp(-0.5*omega.^2).*(abs(omega)<=5);
  Probability(i) = (abs(J(i))/(nasymp*sqrt(2*pi)*model.TruncCorrection))*sum(ProbMC);
end
