function stock_simulation()
% test grn()
%	binner_init(-3,3,30);
%	for i=1:10000
%		binner(grn());
%	end
%	[svals,counts] = binner_counts();
%	hold off;
%	figure(2);
%	plot(svals,counts,'o-');
%stop
		

	s0 = 1
	dt = 1
	T = 250
	mu = 0.12
	sigma = 0.24
	ntrials = 300

	delta_t = dt/T
	nsteps = floor(T/dt)
	data = zeros(nsteps+1,2);
	sfinal = zeros(ntrials,1);
	nbins = 30
	smin = 0
	smax = 2.3
	binner_init(smin,smax,nbins);
	
	hold off;
	figure(1);
	plot([0],[0]);
	hold on;
	for k=1:ntrials
		s = s0;
		data(1,:)=[0,s];
		for i=1:nsteps
%			eps = 2*rand(1)-1;  % uniform random number in [-1,1]

			eps = grn();	% gaussian random number, mean 0, variance 1
			s = s*(1 + mu*delta_t + sigma*eps*sqrt(delta_t));

			data(i+1,:) = [i*dt,s];
		end
		sfinal(k) = s;
		binner(s);
		plot(data(:,1),data(:,2));
		hold on;
	end
	hold off;
%	sfinal
	[svals,counts] = binner_counts();
	figure(2);
	plot(svals,counts,'o-');
	xlabel('final s-value');
	ylabel('frequency');

	[mu,sigma] = binner_stats();
	mu
	sigma

end

function v = grn()
	v = sqrt(2)*erfinv((2*rand()-1));  % pretty slow, it turns out
									   % Compare Box-Mueller!
end

%%%%%%%%%%%% functions below only to make a histogram of sfinal

function binner_init(slo,shi,nbins)
	global myslo myshi mynbins mysrange myds mycounts mycalls mytot mysqtot
	myslo = slo;
	myshi = shi;
	mynbins = nbins;
	mysrange = shi-slo;
	myds = mysrange/nbins;
	mycounts = zeros(nbins,1);
	mycalls = 0;
	mytot = 0;
	mysqtot = 0;
end

function binner(s)
	global myslo myshi mynbins mysrange myds mycounts mycalls mytot mysqtot
	if (s < myslo) || (s > myshi)
		s
		[myslo,myshi]
		fprintf('Warning: binner received data out of bounds\n');
		return
	end
	bin = floor( (s-myslo)*mynbins/mysrange ) + 1;
	mycounts(bin) = mycounts(bin) + 1;
	mycalls = mycalls + 1;
	mytot = mytot + s;
	mysqtot = mysqtot + s^2;
end

function [bincenters,counts] = binner_counts()
	global myslo myshi mynbins mysrange myds mycounts mycalls
	indices = 1:mynbins;
	bincenters = myslo + (indices - 1)*myds + myds/2;
	counts = mycounts; %/mycalls;
end

function [mu,sigma] = binner_stats()
	global mycalls mytot mysqtot
	mu = mytot/mycalls;
	sigma = sqrt( mysqtot/mycalls - mu^2 );
end