sde_library_setup.m

SIMULATION AND ESTIMATION OF STOCHASTIC DIFFERENTIAL EQUATIONS WITH MATLAB

function sde_setup_output = SDE_library_setup(sde_setup_input)

% Library setup: here the number of the state variables, the parameters for
% the models of the library and basic constructs for the library usage are
% defined.
%
% IN:  sde_setup_input; a structure containing the following fields:
%                       sde_setup_input.loaddata; load data from datafile? can be Yes ('Y') or Not ('N')
%                       sde_setup_input.parestimate; estimate parameters? can be Yes ('Y') or Not ('N')
%                       sde_setup_input.time; the array of unique observation times
%                       sde_setup_input.vrbl; the array of unique label-variables
%                       sde_setup_input.xobs; the matrix-shaped observed data
% OUT: sde_setup_output; a structure containing the following fields:
%                        sde_setup_output.bigtheta; complete structural parameter vector
%                        sde_setup_output.parbase; the same as bigtheta, provides parameters starting values for the optimization procedure 
%                        sde_setup_output.parmask; an array containing ones in correspondence of the parameters in bigtheta to be estimated 
%                                                  and zeros in correspondence of the parameters to be held fixed (constant); it has the same 
%                                                  length of bigtheta.                                
%                        sde_setup_output.problem; the name of the current problem/experiment/example etc. (e.g. 'M1')
%                        sde_setup_output.sdetype; the SDE definition: can be 'Ito' or 'Strat' (Stratonovich)
%                        sde_setup_output.numdepvars; the number of dependent variables, i.e. the SDE dimension
%                        sde_setup_output.numsim; the number of desired simulations for the SDE numerical integration 
%                        sde_setup_output.model; the model name (e.g. 'M1a', 'M1b', etc.).
%                        sde_setup_output.owntime; vector containing the equispaced simulation times sorted in ascending order. 
%                                                  It has starting simulation-time in first and ending simulation-time in last position. 
%                                                  Thus OWNTIME(i) - OWNTIME(i-1) = h, where h is the fixed stepsize for the numerical intregration (i=2,3,...)
%                        sde_setup_output.time; the array of unique observation times
%                        sde_setup_output.vrbl; the array of unique label-variables 
%                        sde_setup_output.integrator; the SDE fixed stepsize numerical integration method: can be 'EM' (Euler-Maruyama) or 'Mil' (Milstein)

% Copyright (C) 2007, Umberto Picchini  
% [e-mail]
% http://www.biomatematica.it/Pages/Picchini.html
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
% 
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

MODEL = upper(input('\nWrite the name of the chosen model (e.g. M1a for MODEL 1 with Ito definition, M1b for MODEL 1 with Stratonovich definition, etc.): ','s'));

% For each model the name of the problem (e.g. 'M1a'), the number of
% the state variables (NUMDEPVARS), the SDE type (Ito or Stratonovich) and
% the parameters are specified.
%::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

LOADDATA = sde_setup_input.loaddata;
PARESTIMATE = sde_setup_input.parestimate;
TIME = sde_setup_input.time;
VRBL = sde_setup_input.vrbl;
XOBS = sde_setup_input.xobs;
    
switch MODEL
    
    case 'M1A'
          PROBLEM = 'M1';
          SDETYPE = 'Ito';
          NUMDEPVARS = 1;
          fprintf('\nYou choose     dXt = -a * Xt * dt + sigma * dWt,   X(0) = X0');
          if( (strcmp(LOADDATA,'N') && strcmp(PARESTIMATE,'N')) || (strcmp(LOADDATA,'N') && strcmp(PARESTIMATE,'Y')) ||  (strcmp(LOADDATA,'Y') && strcmp(PARESTIMATE,'N')) )
             a  = input('\n\nWrite the value of the ''a'' parameter: ');
             if(isempty(a))
                error('''a'' must be specified');
             end
             sigma  = input('\nWrite the value of the ''sigma'' parameter: ');
             if(isempty(sigma))
                error('''sigma'' must be specified');
             end
             Xzero = input('\nWrite the value of the initial condition X0: ');
             if(isempty(Xzero))
                error('X0 must be specified');
             end
             % store the parameters into the bigtheta array
             bigtheta(1) = Xzero;  
             bigtheta(2) = a; 
             bigtheta(3) = sigma;
             PARBASE = bigtheta; % the complete array of the user defined parameter values
          elseif( strcmp(LOADDATA,'Y') && strcmp(PARESTIMATE,'Y') )
             fprintf('\n'); 
             bigtheta(1:NUMDEPVARS) = XOBS(1,:);  % subsititute the state variable(s) initial condition(s) with the first observed value(s)
             PARBASE = [];
          end  % if( (strcmp(LOADDATA,'N') || strcmp(LOADDATA,'n')) && (strcmp(PARESTIMATE,'Y') || strcmp(PARESTIMATE,'y')) )
          %:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
          PARMASK = [0,1,1];  % write 1 for parameters to be estimated and 0 for fixed parameters. WARNING: PARMASK(1) should always be set to zero (corresponds to the SDE initial condition)
          %:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    case 'M1B' 
          PROBLEM = 'M1';
          SDETYPE = 'Strat';
          NUMDEPVARS = 1;
          fprintf('\n\nYou choose     dXt = -a * Xt * dt + sigma o dWt,   X(0) = X0');
          if( (strcmp(LOADDATA,'N') && strcmp(PARESTIMATE,'N')) || (strcmp(LOADDATA,'N') && strcmp(PARESTIMATE,'Y')) ||  (strcmp(LOADDATA,'Y') && strcmp(PARESTIMATE,'N')) )
             a  = input('\n\nWrite the value of the ''a'' parameter: ');
             if(isempty(a))
                error('''a'' must be specified');
             end
             sigma  = input('\nWrite the value of the ''sigma'' parameter: ');
             if(isempty(sigma))
                error('''sigma'' must be specified');
             end
             Xzero = input('\nWrite the value of the initial condition X0: ');
             if(isempty(Xzero))
                error('X0 must be specified');
             end
             % store the parameters into the bigtheta array
             bigtheta(1) = Xzero;  
             bigtheta(2) = a; 
             bigtheta(3) = sigma;     
             PARBASE = bigtheta; % the complete array of the user defined parameter values
          elseif( strcmp(LOADDATA,'Y') && strcmp(PARESTIMATE,'Y') )
             fprintf('\n');
             bigtheta(1:NUMDEPVARS) = XOBS(1,:);  % subsititute the state variable(s) initial condition(s) with the first observed value(s)
             PARBASE = [];
          end  % if( (strcmp(LOADDATA,'N') || strcmp(LOADDATA,'n')) && (strcmp(PARESTIMATE,'Y') || strcmp(PARESTIMATE,'y')) )
          %:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
          PARMASK = [0,1,1];  % write 1 for parameters to be estimated and 0 for fixed parameters. WARNING: PARMASK(1) should always be set to zero (corresponds to the SDE initial condition)
          %:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
          
     %-------------------------------------------------------------------------------------------------------------------
     
     case 'M2A' 
          PROBLEM = 'M2';
          SDETYPE = 'Ito';
          NUMDEPVARS = 1;
          fprintf('\n\nYou choose     dXt = (a * Xt + b) * dt + sigma * dWt,   X(0) = X0');
          if( (strcmp(LOADDATA,'N') && strcmp(PARESTIMATE,'N')) || (strcmp(LOADDATA,'N') && strcmp(PARESTIMATE,'Y')) ||  (strcmp(LOADDATA,'Y') && strcmp(PARESTIMATE,'N')) )
             a  = input('\n\nWrite the value of the ''a'' parameter: ');
             if(isempty(a))
                error('''a'' must be specified');
             end
             b  = input('\nWrite the value of the ''b'' parameter: ');
             if(isempty(b))
                error('''b'' must be specified');
             end
             sigma  = input('\nWrite the value of the ''sigma'' parameter: ');
             if(isempty(sigma))
                error('''sigma'' must be specified');
             end
             Xzero = input('\nWrite the value of the initial condition X0: ');
             if(isempty(Xzero))
                error('X0 must be specified');
             end
            % store the parameters into the bigtheta array
             bigtheta(1) = Xzero;  
             bigtheta(2) = a; 
             bigtheta(3) = b;
             bigtheta(4) = sigma;
             PARBASE = bigtheta; % the complete array of the user defined parameter values
          elseif( strcmp(LOADDATA,'Y') && strcmp(PARESTIMATE,'Y') )
             fprintf('\n');
             bigtheta(1:NUMDEPVARS) = XOBS(1,:);  % subsititute the state variable(s) initial condition(s) with the first observed value(s)
             PARBASE = [];
          end  % if( (strcmp(LOADDATA,'N') || strcmp(LOADDATA,'n')) && (strcmp(PARESTIMATE,'Y') || strcmp(PARESTIMATE,'y')) )
          %:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
          PARMASK = [0,1,1,1];  % write 1 for parameters to be estimated and 0 for fixed parameters. WARNING: PARMASK(1) should always be set to zero (corresponds to the SDE initial condition)
          %:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
     case 'M2B' 
          PROBLEM = 'M2';
          SDETYPE = 'Strat';
          NUMDEPVARS = 1;
          fprintf('\n\nYou choose     dXt = (a * Xt + b) * dt + sigma o dWt,   X(0) = X0');
          if( (strcmp(LOADDATA,'N') && strcmp(PARESTIMATE,'N')) || (strcmp(LOADDATA,'N') && strcmp(PARESTIMATE,'Y')) ||  (strcmp(LOADDATA,'Y') && strcmp(PARESTIMATE,'N')) )
             a  = input('\n\nWrite the value of the ''a'' parameter: ');
             if(isempty(a))
                error('''a'' must be specified');
             end
             b  = input('\nWrite the value of the ''b'' parameter: ');
             if(isempty(b))
                error('''b'' must be specified');
             end
             sigma  = input('\nWrite the value of the ''sigma'' parameter: ');
             if(isempty(sigma))
                error('''sigma'' must be specified');
             end
             Xzero = input('\nWrite the value of the initial condition X0: ');
             if(isempty(Xzero))
                error('X0 must be specified');
             end
            % store the parameters into the bigtheta array
             bigtheta(1) = Xzero;  
             bigtheta(2) = a; 
             bigtheta(3) = b;
             bigtheta(4) = sigma;   
             PARBASE = bigtheta; % the complete array of the user defined parameter values
          elseif( strcmp(LOADDATA,'Y') && strcmp(PARESTIMATE,'Y') )
             fprintf('\n');
             bigtheta(1:NUMDEPVARS) = XOBS(1,:);  % subsititute the state variable(s) initial condition(s) with the first observed value(s)
             PARBASE = [];
          end  % if( (strcmp(LOADDATA,'N') || strcmp(LOADDATA,'n')) && (strcmp(PARESTIMATE,'Y') || strcmp(PARESTIMATE,'y')) )
          %:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
          PARMASK = [0,1,1,1];  % write 1 for parameters to be estimated and 0 for fixed parameters. WARNING: PARMASK(1) should always be set to zero (corresponds to the SDE initial condition)
          %:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
          
     %-------------------------------------------------------------------------------------------------------------------
     
     case 'M3A' 
          PROBLEM = 'M3';
          SDETYPE = 'Ito';
          NUMDEPVARS = 1;
          fprintf('\n\nYou choose     dXt = (a - sigma^2/2) * dt + sigma * dWt,   X(0) = X0');
          if( (strcmp(LOADDATA,'N') && strcmp(PARESTIMATE,'N')) || (strcmp(LOADDATA,'N') && strcmp(PARESTIMATE,'Y')) ||  (strcmp(LOADDATA,'Y') && strcmp(PARESTIMATE,'N')) )
             a  = input('\n\nWrite the value of the ''a'' parameter: ');
             if(isempty(a))
                error('''a'' must be specified');
             end
             sigma  = input('\nWrite the value of the ''sigma'' parameter: ');
             if(isempty(sigma))