function [model,interp,theta,trendparams] = InitLogLikPb(Pobs,params,r,delta,mina,maxa,N,trend,options)
% INITLOGLIGPB Initialization of the problem to be maximized

%% Initialization
T         = length(Pobs);
t         = (1:1:T)'/T;
switch trend
  case 'without'
    lengthfg = 0;
  case 'linear'
    lengthfg = 1;
  case {'quadratic','RCS3'}
    lengthfg = 2;
  case {'cubic','RCS4'}
    lengthfg = 3;
end
if isempty(params), params = {'first-guess' [zeros(1,lengthfg) 1 0.3 0.02]}; end

%% First guess
switch params{1}
  case 'first-guess'
    fgparams     = params{2};
    switch trend
      case 'without'
        g            = zeros(T,1);
        trendparams  = zeros(0,1);
      case {'linear','quadratic','cubic'}
        poly           = Legendre(t);
        lPobslm        = fitlm(poly(:,1:lengthfg),log(Pobs));
        trendparams    = lPobslm.Coefficients.Estimate(2:end)+...
                         lPobslm.Coefficients.SE(2:end).*fgparams(1:lengthfg)'; % fgparams: normal distrib
        g              = poly(:,1:lengthfg)*trendparams;
        model.basis    = poly(:,1:lengthfg);
      case {'RCS3','RCS4'}
        model.basis    = importdata(fullfile('Basis',['PricensBSM' num2str(lengthfg+1) '.mat']));
        lPobslm        = fitlm(model.basis,log(Pobs));
        trendparams    = lPobslm.Coefficients.Estimate(2:end)+...
                         lPobslm.Coefficients.SE(2:end).*fgparams(1:lengthfg)'; % fgparams: normal distrib
        g              = model.basis*trendparams;
    end
    Psto         = Pobs.*exp(-g);
    EP           = mean(Psto);
    sigmaP       = std(Psto);
    Prior        = [fgparams(end-2)*EP -sigmaP/fgparams(end-1) delta fgparams(end)*EP]';
    model.params = [Prior(1:4); r];
    theta        = [trendparams; Prior];
  case 'params'
    theta        = params{2};
    model.params = [theta((end-3):end); r];
    switch trend
      case {'without','linear','quadratic','cubic'}
        poly         = Legendre(t);
        model.basis  = poly(:,1:lengthfg);
      case {'RCS3','RCS4'}
        model.basis  = importdata(fullfile('Basis',['PricensBSM' num2str(lengthfg+1) '.mat']));
    end
end

%% Define the shocks structure by a discretization of the gaussian density
truncated_normal_rule(3,10,0,1,-5,5,fullfile('..','Results','GaussianQuadrature'))
model.shocks.e  = load(fullfile('..','Results','GaussianQuadrature_x.txt'));
model.shocks.w  = load(fullfile('..','Results','GaussianQuadrature_w.txt'));
TruncatedNormal = truncate(makedist('Normal'),-5,5);
model.funrand   = @(nrep,nper) TruncatedNormal.random(nrep,nper);
nrep            = 10000;
nper            = 120;
model.harvest   = [NaN(nrep,1) model.funrand(nrep,nper-1)];
model.TruncCorrection = 2*normcdf(5)-1;

%% Define the grid of interpolation (equally space points)
interp.s        = linspace(mina,maxa,N)';

%% Solve for the REE
interp        = SolveStorageDL(model,interp,options);
