function [model,interp,Psample,target,Inittheta,Nsample,options,N,mina,maxa,r,delta] = InitMonteCarlo(target,type,quartic)

%% Define the grid of interpolation (equally space points)
N                  = 1000; % number of points on the approximation grid
mina               = -5;   % lower limit of the grid
maxa               = 20;   % upper limit of the grid
r                  = 0.05; % Interest rate assumed exogenous and fixed
delta              = 0;    % Decay rate assumed exogenous and fixed at zero (Cafiero and Wright)
interp.s           = linspace(mina,maxa,N)';
interp.PriceInterp = [];

%% 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);
model.TruncCorrection = 2*normcdf(5)-1;

%% Options
options   = struct('ActiveParams', [1 1 0 1],...
                   'bounds'      , struct('lb',[-inf(3,1); 0],...
                                          'ub',[+inf; -1E-3; +inf(2,1)]),...
                   'cov'         , 2,...
                   'Display'     , 0,...
                   'estimator'   , 'UML',...
                   'InterpMethod', 'spline',...
                   'MaxIter'     , 1E3,...
                   'ncpu'        , 8,...
                   'numhessianoptions',struct('FinDiffRelStep'  , 1E-3,...
                                              'UseParallel'     , true),...
                   'seed'        , 1,...
                   'TolX'        , 1E-10,...
                   'UseParallel' , 'never',...
                   'Vectorized'  ,'off');

fminconoptions   = optimoptions('fmincon'      , ...
                                'Algorithm'    , 'sqp',...
                                'DiffMinChange', eps^(1/3),...
                                'Display'      , 'off',...
                                'FinDiffType'  , 'central',...
                                'MaxFunEvals'  , 1E4,...
                                'TolFun'       , 1e-6,...
                                'TolX'         , 1e-6,...
                                'UseParallel'  , true);

fminsearchoptions = optimset('Display'      , 'off',...
                             'MaxFunEvals'  , 1E4,...
                             'MaxIter'      , 1E3,...
                             'TolFun'       , 1e-6,...
                             'TolX'         , 1e-6);

patternsearchoptions = psoptimset('Display'       , 'off',...
                                  'MaxFunEvals'   , 1E5,...
                                  'MaxIter'       , 1E4,...
                                  'CompletePoll'  , 'on',...
                                  'CompleteSearch', 'on',...
                                  'SearchMethod'  , 'searchneldermead',...
                                  'TolFun'        , 1e-6,...
                                  'TolX'          , 1e-6,...
                                  'UseParallel'   , false);

fminuncoptions    = optimoptions('fminunc'      ,...
                                 'Algorithm'    ,'quasi-newton',...
                                 'DiffMinChange', eps^(1/3),...
                                 'Display'      , 'off',...
                                 'FinDiffType'  , 'central',...
                                 'MaxFunEvals'  , 1E4,...
                                 'TolFun'       , 1e-6,...
                                 'TolX'         , 1e-6);

options.solver        = {'fminsearch' 'patternsearch' 'fmincon' 'fminunc'};
options.solveroptions = {fminsearchoptions patternsearchoptions fminconoptions fminuncoptions};

%% Simulation of artificial price samples
rng(options.seed);
Nsample         = 3000;           % Number of replications in the MC experiment
burninper       = 50;             % Number of burn-in periods

if nargin<2 || isempty(type), type = 'long'; end
switch type
  case 'short'
    samplelen = 50;
  case 'long'
    samplelen = 100;
end
nper             = samplelen+burninper;
A0               = zeros(Nsample,1);     % Initial vector of availability
model.params     = [target(1:4) r];
model.basis      = zeros(samplelen,0);

if nargin ==3 ,options.InterpMethod = quartic; end
switch lower(options.InterpMethod)
    case 'quartic'
       [interp,~,pstar] = SolveStorageQuartic(model,interp,options);
       [~,Xsim]         = SimulStorageQuartic(model,interp,A0,nper);
    otherwise
       [interp,~,pstar] = SolveStorageDL(model,interp,options);
       [~,Xsim]         = SimulStorage(model,interp,A0,nper);
end

Psample          = squeeze(Xsim(:,2,burninper+1:end));
target           = [target pstar];

%% Random starting values between -20% and +20% of target values
rng(options.seed);
Inittheta      = zeros(Nsample,5);
Inittheta(:,1) = target(1)*(0.8+0.4*rand(Nsample,1));
Inittheta(:,2) = target(2)*(1.2-0.4*rand(Nsample,1));
Inittheta(:,3) = delta;
Inittheta(:,4) = target(4)*(0.8+0.4*rand(Nsample,1));
Inittheta(:,5) = r;

%% Random supply shock to calculate the marginal distribution
rng(options.seed);
nrep             = 1E4;
nburn            = 20;
nper             = 100+nburn;
harvest          = NaN(nrep,nper);
harvest(:,2:end) = model.funrand(nrep,nper-1);
model.harvest    = harvest;
