sde_library_setup.m

来自「SIMULATION AND ESTIMATION OF STOCHASTIC 」· M 代码 · 共 853 行 · 第 1/4 页

            % store the parameters into the bigtheta array
             bigtheta(1) = Xzero1;   
             bigtheta(2) = Xzero2;  
             bigtheta(3) = alpha1; 
             bigtheta(4) = alpha2; 
             bigtheta(5) = beta11; 
             bigtheta(6) = beta12; 
             bigtheta(7) = beta21; 
             bigtheta(8) = beta22;
             bigtheta(9) = sigma1;
             bigtheta(10)= sigma2;  
             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,0,1,1,1,1,1,1,1,1];  % write 1 for parameters to be estimated and 0 for fixed parameters. WARNING: PARMASK(1) and PARMASK(2) should always be set to zero (correspond to the SDE initial conditions)
          %:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
          
     %-------------------------------------------------------------------------------------------------------------------
     
     case 'M10A'
          PROBLEM = 'M10';
          SDETYPE = 'Ito';
          NUMDEPVARS = 1;
          fprintf('\nYou choose     dXt = a * (b - Xt) * dt + sigma * sqrt(Xt) * 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 'M10B' 
          PROBLEM = 'M10';
          SDETYPE = 'Strat';
          NUMDEPVARS = 1;
          fprintf('\n\nYou choose     dXt = [a * (b - Xt) -1/4 * sigma^2] * dt + sigma * sqrt(Xt) 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)
          %:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
          
     %-------------------------------------------------------------------------------------------------------------------
     
    otherwise
          error('The model name ''%s'' is not included in the list of available models...please check.',MODEL);
end
%:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::








%::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
%                                  DO NOT MODIFY THE CODE BELOW   *    DO NOT MODIFY THE CODE BELOW   *    DO NOT MODIFY THE CODE BELOW    *    DO NOT MODIFY THE CODE BELOW
%::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
%::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
%                                  DO NOT MODIFY THE CODE BELOW   *    DO NOT MODIFY THE CODE BELOW   *    DO NOT MODIFY THE CODE BELOW    *    DO NOT MODIFY THE CODE BELOW
%::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
%::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
%                                  DO NOT MODIFY THE CODE BELOW   *    DO NOT MODIFY THE CODE BELOW   *    DO NOT MODIFY THE CODE BELOW    *    DO NOT MODIFY THE CODE BELOW
%::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
%::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
%                                  DO NOT MODIFY THE CODE BELOW   *    DO NOT MODIFY THE CODE BELOW   *    DO NOT MODIFY THE CODE BELOW    *    DO NOT MODIFY THE CODE BELOW
%::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
%::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
%                                  DO NOT MODIFY THE CODE BELOW   *    DO NOT MODIFY THE CODE BELOW   *    DO NOT MODIFY THE CODE BELOW    *    DO NOT MODIFY THE CODE BELOW
%::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::





% General settings
NUMSIM = input('\nNumber of trajectories to be simulated: ');
[n1,n2]=rat(NUMSIM);
if(~isequal(n2,1) || isempty(NUMSIM) || NUMSIM <= 0)
    error('The number of trajectories NUMSIM must be a positive integer');
end

if(strcmp(LOADDATA,'N'))
   T0 = input('\nStarting time T0: ');
   if(isempty(T0))
       error('T0 must be specified');
   end
   T  = input('\nFinal time T: ');
   if(isempty(T))
       error('T must be specified');
   end
   if(T<=T0)
       error('T must be greater than T0')
   end
   fprintf('\nChoose the integration stepsize (should be << T - T0):');
   h  = input('\nh: ');
   if(h<=0 || isempty(h))
       error('The stepsize must be a positive number') 
   end 
   if(h>(T-T0))
       error('Choose a smaller stepsize')
   end
elseif (strcmp(LOADDATA,'Y'))
       T0 = min(TIME);
       T = max(TIME);
       fprintf('\nChoose the integration stepsize (should be << %d):', min(diff(TIME(diff(TIME)>0))));
       h  = input('\nh: ');
       if(h > min(diff(TIME(diff(TIME)>0))))
          error('Choose a smaller stepsize')
       end
       if(h<=0 || isempty(h))
          error('The stepsize must be a positive number') 
       end    
end


%:::::: the desired discretization for [T0,T] :::::::::::::::::::::
OWNTIME = T0:h:T;    
if(OWNTIME(end)~=T)
   OWNTIME(end)=T;
end
%::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

if(strcmp(SDETYPE,'Ito'))
   fprintf('\nChoose an integration method:');
   fprintf('\nMethods for Ito SDE: Euler-Maruyama (EM), Milstein (Mil);');
   INTEGRATOR = upper(input('\n''EM'' or ''Mil''?: ','s'));
   if(strcmp(INTEGRATOR,'EM')==0 && strcmp(INTEGRATOR,'MIL')==0)
       error('You must specify ''EM'' or ''Mil''.');
   end
else
   INTEGRATOR = 'MIL';
end

if (strcmp(LOADDATA,'N') && strcmp(PARESTIMATE,'Y'))
    NTIMEOBS  = input('\nNumber of observation time-points: '); 
    [n1,n2]=rat(NTIMEOBS);
    if(~isequal(n2,1) || isempty(NTIMEOBS) || NTIMEOBS <=0)
       error('The number of observation time-points must be a positive integer');
    end
    delta = (T-T0)/(NTIMEOBS-1);
    TIME = T0:delta:T;   % create NTIMEOBS equally spaced observation-times into [T0,T] (TIME should be such that length(TIME)<<length(OWNTIME))
    if(TIME(end)~=T)
       TIME(end)=T;
    end
    % Now replicate the values in TIME depending on the NUMDEPVARS value (useful for future versions of SDE_Toolbox)
    TIME = TIME(ones(1,NUMDEPVARS),:);
    TIME = TIME(:).';
    VRBL = [1:NUMDEPVARS];
    % Now replicate the values in VRBL depending on the TIME value (useful for future versions of SDE_Toolbox)
    VRBL = repmat(VRBL,1,length(unique(TIME)));
end
                     

sde_setup_output = struct('bigtheta',bigtheta,'parbase',PARBASE,'parmask',PARMASK,'problem',PROBLEM,'sdetype',SDETYPE,'numdepvars',NUMDEPVARS,'numsim',NUMSIM,'model',MODEL,'owntime',OWNTIME,'time',TIME,'vrbl',VRBL,'integrator',INTEGRATOR);