function batchruns(batchnum,noload)
  % batchnum is a numeric ID to identify a given model specification
  % if noload = 1, then do NOT load config.  Instead start fresh and overwrite the config file for this batchnum.
  
  if nargin==0;  batchnum=1; noload = 0;  end;  % can either pass in batchnum or hardcode it
  if nargin==1;  noload = 0;  end;

  noload = 1;  % don't load existing config file
  pause off;

% Set which batch job to start with. If the relevant config file doesn't
% exist, use the default values and save to the config file.

% Build Par struct with general program data
Par.batchjob        = batchnum;
Par.small           = 1e-8;
Par.nocheck         = -1e50; 
Par.uu_competition  = 1;
Par.uu_monopoly     = 2;
Par.uu_planner      = 3;
Par.uu_quality      = 1;
Par.uu_cost         = 2;
Par.uu_capacity     = 3;
Par.uu_ran_entry    = 1;
Par.uu_det_entry    = 2;


deltagrid = ( 0:.1:.9 )';      % construct comparative static in delta, as an example
ngrid = length(deltagrid);
innov = zeros(ngrid,1);

for job = 1:ngrid;

  Par.batchjob = Par.batchjob + 1;
  Par.configfile = ['config_' int2str(Par.batchjob) '.mat'];

  % Load Params struct from file (if it exists) or use the default Params
  if noload==0 && exist(Par.configfile, 'file')
    disp('loading config file...');
    load(Par.configfile, 'Params');
  else
    disp('starting with default config file...');
    Params = default(Par);
    save(Par.configfile, 'Params');
  end
  
  Params.delta = deltagrid(job);        %% hard code parameters here if don't want values set below in  default() 
  Params.rlg_foc = 1;
  Params.kmax = 9;
  Params.entry_at = -3;
  Params.max_firms = 3;
  Params.rlg_outgood = 1;  Params.rlg_y = 15; Params.rlg_inv = .2;  % Params.rlg_leap = .2;
  
  % Params.wstar = 7;  Params.rlg_y = 0;  Params.rlg_outgood = 1;
  
  save(Par.configfile, 'Params');
  
  
  % Display current Params
  showparms(Par, Params);
  
  
  profit_main(Par.batchjob);  % Get profits
  
  eql_ma_foc(Par.batchjob);   %  Get Equilibrium
  
  simstuff = welf_ma(Par.batchjob);      % Simulate model
  
  innov(job,1) = mean( simstuff.totinnov );

end;  %% end for loop constructing a comparative static

disp(' ');
disp('------  Comparative Static for Innovation as function of delta  --------');
disp(' ');
disp('       delta   industry innovation');
fprintf(' %12.4f  %12.4f \n', [deltagrid innov]');

%  HEY:  Should probably modify  eql_ma_foc() and welf_ma() to return the items of interest for the comparative static.  
%  Would be easiest to return in a structure.

end % End of function batchruns



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function default - load default values into Param struct
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function p = default(Par)
p.eql_type     = Par.uu_competition;    % = Par.uu_competition  OR  Par.uu_monopoly  OR  Par.uu_planner, which are constants: 1, 2, 3, respectively
p.ind_type     = Par.uu_quality;        % = Par.uu_quality  OR  Par.uu_cost  OR  Par.uu_capacity,  which are constants: 1, 2, 3, respectively
p.entry_type   = Par.uu_ran_entry;      % = Par.uu_ran_entry  OR  Par.uu_det_entry,  which are constants: 1, 2, respectively
p.max_firms    = 4;
p.entry_low    = 22;                    % lower bound of uniform entry cost distribution
p.entry_high   = 24;                    % upper bound of uniform entry cost distribution
p.entry_sunk   = 0.2;                   % fixed (non-random) entry cost (not recommended)
p.scrap_val    = 21;                    % lower bound of scrap value
p.rlg_scrap_hi = 22;                    % upper bound of scrap value
p.entry_at     = 4;                     % efficiency level of entrant (1 to kmax)
p.beta         = 0.925;                 % discount factor
p.delta        = 0.5;                   % probability outside good improves
p.inv_mult     = .5;                    % investment efficiency parameter a in  ax/(1+ax)  for prob of success
p.mc           = 5;                     % constant marginal cost of production
p.mkt_size     = 5;                     % market size
p.kmax         = 7;                     % max firm efficiency (i.e., # rungs on quality ladder)
p.wstar        = -1;                    % if > 0, use PM model with wstar as point of curvature that bounds state space by reducing incentives for leaders to innovate
                                        % if wstar = -1, use linear model of GG
					
p.intercept    = 3;                     % demand intercept.  Only relevant for COST or CAPACITY investment game, not Quality Ladder game
p.fixed_cost   = 0.0;                   % fixed cost per period
p.gamma        = 1;                     % maximum marginal cost in COST investment game
p.tau          = 0.1;                   % capacity ratio in CAPACITY investment game

p.rlg_wscale   = .2;                    % quality coefficient (3 in RAND 1994)
p.rlg_wshift   = 0;                     % shift in firms' (scaled) qualities, relative to outside good at zero (-7 in RAND 1994)
p.rlg_sh_cap   = 1;                     % maximum share of lead firm:  .55 and .65 counterfactuals in RAND 1994
p.rlg_inv      = 0;                     % spillover: if >0, all laggards have investment multiplier of  inv_mult*(1+spillover).  If <0, then inv_mult*(1 + |spillover| * #Steps_Behind)
p.rlg_y        = 15;                    % Income: if 0, then linear utility.  >0 for Cobb-Douglas utility = w+log(y-p)-log(y)  which bounds p<y
p.rlg_alpha    = .4;                    % price coefficient
p.rlg_outgood  = 1;                     % if 1, outside good present.  if 0, no outside good
p.rlg_leap     = 0.0;                   % Prob entrant goes to KMAX instead of ENTRY_AT.  Leader who successfully innovates, will have a 1-step lead over a KMAX entrant since entry occurs in model before investment outcome
p.rlg_maxp     = 9999;                  % ignored unless NO Outside Good  &  Y = 0, in which case monopolist price = rlg_maxp * price of leader in max-differentiated duopoly
p.rlg_min_innov1 = 0;                   % minimum innovation rate by frontier firms.  Should generally be 0.  No justification for >0.

p.rlg_foc      = 0;                     % if 1, use FOC method.  If 0, use Iterative Best Response method of Pakes-McGuire

p.profit_done  = 0;                     % flag for profit matrix computed
p.eql_done     = 0;                     % flag for equilibrium computed
p.active_cfg   = Par.batchjob;          % batchjob number.  Not sure why we have both  Par.batchjob and Params.active_cfg

end % End of function default



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function acronym - get textual description of types
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [s,  etype,  itype] = acronym(Par, et, it)
if et == Par.uu_competition
    s = '';
    etype = 'markov perfect nash equilibrium';
elseif et == Par.uu_monopoly
    s = 'm';
    etype = 'monopoly';
elseif et == Par.uu_planner
    s = 's';
    etype = 'social planner';
end

if it == Par.uu_quality
    s = [s 'b'];
    itype = 'quality (differentiated products)';
elseif it == Par.uu_cost
    s = [s 'c'];
    itype = 'marginal cost (homogenous products)';
elseif it == Par.uu_capacity
    s = [s 'p'];
    itype = 'capacity (homogenous products)';
end

if Par.batchjob > 0 
    s = ['bj' int2str(Par.batchjob) '_' s]; 
end
end % End of function acronym



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function showparms - Display current parameters in menu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function showparms(Par, Params)

[~,et,it] = acronym(Par,Params.eql_type, Params.ind_type);

disp('********                          Dynamic Oligopoly Algorithm                          ********');
disp('******** based on Goettler and Gordon (2013) modification of Pakes and McGuire (1994)  ********');
cfg = ['Current configuration: ' Params.active_cfg];
disp( center(cfg, ''));

disp( outline('Model Primitives'));
fprintf(['INVESTMENT IN:\t' it '\n']);
fprintf(['EQUILIBRIUM:\t' et '\n']);
fprintf(['MAX NUMBER OF FIRMS:\t' int2str(Params.max_firms) '\tMAX EFFICIENCY LEVEL (KMAX):\t' int2str(Params.kmax) '\n']);
disp( outline('Static Parameters'));
if (Params.ind_type == Par.uu_quality)
    fprintf(['MARGINAL, FIXED COST:\t' num2str(Params.mc) '\t' num2str(Params.fixed_cost) '\t']);
    fprintf(['MKT SIZE:\t' num2str(Params.mkt_size) '\n']);
    fprintf(['OUTSIDE GOOD:\t' int2str(Params.rlg_outgood) '\t']);
    fprintf(['Y:\t' num2str(Params.rlg_y) '\t']);
    fprintf(['MAXP FACTOR: ' num2str(Params.rlg_maxp) '\n']);
    fprintf(['W SCALE,  SHIFT:\t' num2str(Params.rlg_wscale) '\t' num2str(Params.rlg_wshift) '\t']);
    fprintf(['PRICE COEFF:\t' num2str(Params.rlg_alpha) '\t']);
    fprintf(['SHARE CAP:  ' num2str(Params.rlg_sh_cap) '\n']);
    fprintf(['w*:\t' int2str(Params.wstar) '  (-1 for linear GG, >0 PM w*)\n']);
    if Params.rlg_outgood==0;
      if Params.rlg_y>0;  disp('Note: with no outside good & Y>0 , monopoly price is Y');  
      else;               disp('Note: with no outside good & Y=0 , monopoly price is MAXP * leader price in max-differentiated duopoly');
      end;
    end;
    % if Params.wstar>0 ;
    %  disp(['Note: w*= ' num2str(Params.wstar) ' yields log(eg(kmax)) = ' num2str(log(exp(Params.wstar)*(2.0-exp(-(Params.kmax*Params.rlg_wscale+Params.rlg_wshift-Params.wstar)))))]);
    % end;
elseif (Params.ind_type == Par.uu_cost)
    disp(['demand intercept: ' num2str(Params.intercept)]);
    disp('');
    disp(['fixed cost: ' num2str(Params.fixed_cost)]);
    disp('');
    disp(['maximum marginal cost: ' num2str(Params.gamma)]);
elseif (Params.ind_type == Par.uu_capacity)
    disp(['demand intercept: ' num2str(Params.intercept)]);
    disp('');
    disp(['marginal cost: ' num2str(Params.mc)]);
    disp('');
    disp(['tau: ' num2str(Params.tau)]);
end

disp( outline('Dynamic Parameters'));
fprintf(['DISCOUNT FACTOR:\t' num2str(Params.beta) '\tPROB OUTSIDE GOOD IMPROVES (delta): ' num2str(Params.delta) '\n']);
fprintf(['SCRAP VALUE: ' num2str(Params.scrap_val) '-' num2str(Params.rlg_scrap_hi) '\t']);
fprintf('SUNK COST: ');
if Params.entry_type == Par.uu_det_entry;
  fprintf(['fixed at ' num2str(Params.entry_sunk) '\t']);
else;
  fprintf([num2str(Params.entry_low) '-' num2str(Params.entry_high) '\t']);
end;
fprintf(['ENTRY AT:' int2str(Params.entry_at) ' (if <0, relative to leader) \n']);
fprintf(['INNOVATION EFFICIENCY: a = ' num2str(Params.inv_mult) '  spill = ' num2str(Params.rlg_inv) '\tPROB ENTER AT KMAX: ' num2str(Params.rlg_leap) '\n']);
disp( repmat('-', 1,79));
end % End of function showparms



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Helper functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% centered formatting - for headings
function sfmt = center(s,  fill)
spad = repmat(fill, 1, floor( (80-(length(s)+2)) /2) );
sfmt = [ spad  ' ' s ' ' spad];
end % End of function center

% outlined formatting - for minor headings
function sfmt = outline(title)
sfmt = ['---( ', title, ' )'];
sfmt = [sfmt repmat('-', 1,79-length(sfmt))];
end % End of function outline


% get a number from the user
function ret = getnum(prompt,  default)
ret = input(prompt);
if (length(ret) == 0)
    ret = default;
end
end % End of function getnum
