Feeds:
Posts

## REMINISCES OF LEHMAN BROTHERS FIXED INCOME RESEARCH 1995

My first job was at Lehman Brothers Fixed Income Research in Andy Morton’s group.  I was hard-working and naive about what was required then for vertical ascent.  I came with mathematical background but untrained in science and I spent my time coding.  My perspective has matured regarding science since, so now I realize the hard work of getting nature to follow the quantitative models now attempting to get a static Einstein universe to produce the CMB anisotropy.  What is clear is that the hard part of numerical work remains as much drudgery as I had in the trading floor of Lehman, 13 hours in the neon light punctuated by dark skies.  What has not changed when making Octave models for CMB anisotropy (as opposed to the derivative pricing models) is that the annoyance of the coding problems is dwarfed by the difficulty of producing a model that can fit the observed anisotropy.  Tonight, I can’t solve the problem: find a distribution of point masses on a three-sphere such that the gravitational redshifts alone may reproduce the anisotropy.  So much for the coding problem.

The sample-frequency-concentrated-distribution looks quite good but not in a realistic FREQUENCY SCALE.

```pkg load statistics

global n=10000
global q0=100
global points=zeros(n,2)
global q0
global n
global X = mvnrnd( [0,0,0,0],diag([1,1,1,1]),n);

for i=1:n,
X(i,:)=X(i,:)/norm(X(i,:),2);
endfor

# stereographic projection + normalize

function z=pwrrnd(k,n)
u=unifrnd(0,1,n)
z=(1-u)^(1/(1-k))*0.001
endfunction

global masses = lognrnd(0,4,n,1)*50;

function d=direction(x1,x2)
d=x2-x1;
d=d/norm(d,2);
endfunction

function d=distance(x1,x2)
d=acos(dot(x1,x2));
endfunction

function f=gforce( x1,x2, m1,m2)
f=m1*m2*direction(x1,x2)/distance(x1,x2)^2;
#f=0
endfunction

timestep=10

if 0,
for a=1:n,
force_on_x = 0;

x1 = X(a,:);
m1 = masses(a);
for b=1:n,
x2 = X(b,:);
d=distance(x1,x2);
m2=masses(b);
if d>0 && d<pi,
force_on_x = force_on_x + gforce(x1,x2,m1,m2);
endif
endfor
x = x+0.5*force_on_x*timestep^2/m1;
x = x/norm(x,2);
X(a,:)=x;
endfor
endif

global z = zeros(n,3);
for i=1:n,
y=X(i,:);
y1=y(1);
y2=y(2);
y3=y(3);
y4=y(4);
zp=[ y1/(1-y4), y2/(1-y4),y3/(1-y4) ];
z(i,:) = zp/norm(zp,2);

# select from log-normal
#u = lognrnd(0,1);
#if u>5,
#z(i,:) = [0,0,0];
#endif
endfor

function t=fitangle(x)
t=x;
while t>2*pi,
t=t-2*pi;
endwhile
while t<0,
t=t+2*pi;
endwhile
endfunction

function b=bin_s2point( z )
if abs(z)<0.001,
b=[0,0];
return
endif
global q0;
theta = atan( z(2)/z(1));
theta=fitangle(theta);
phi = atan( sqrt(z(1)^2+z(2)^2)/z(3));
phi = fitangle(phi);
thetan = floor(theta*q0/(2*pi));
phin = floor(phi*q0/(2*pi));
b=[thetan,phin];
endfunction
# create a spherical grid in longitudes/latitudes
global g = zeros(q0,q0);

global refpoint = [1,0,0,0]

function dox()
global n;
global g;
global X;
global masses;
for i = 1:n,
global refpoint;
global z;
global q0;
global points;
dg = zeros(q0,q0);

cds = bin_s2point(z(i,:));
cds(1)=mod(cds(1),q0);
cds(2)=mod(cds(2),q0);
disp([cds(1),cds(2)])
points(i,:)=cds;
fflush(stdout)
#sdir = sign(direction(refpoint,X(i,:));
dist = distance(refpoint,X(i,:));
#g(cds(1)+1,cds(2)+1) = g(cds(1)+1,cds(2)+1)+1;
m1=masses(i);
dg(cds(1)+1,cds(2)+1) = m1/dist;
dg = imsmooth( dg, "Gaussian", ceil(2),ceil(2));
#disp(sum(sum(dg)))
#fflush(stdout)
dg(1,1)=0;
g = g+dg;
endfor
endfunction

dox()
sum(sum(g))
#g = (g-min(min(g)));
g=g/sum(sum(g));

# value at a spherical harmonic

function v=sphCoeff( l,m, g)
v = 0;
q0=20;
for j=1:q0,
for k=1:q0,
w1 = cos(2*pi*j/q0);
w2 = exp(1i*2*pi*k/q0);
fnval = gsl_sf_legendre_sphPlm( l, m, w1 )*w2;
if ~isnan(fnval*g(j,k)),
v = v + fnval*g(j,k);
endif
endfor
endfor
disp([l,m,v])
fflush(stdout)
endfunction

function v=sumSphCoeff( l, g)
v=0;
for m=0:1
v = v+sphCoeff(l,m,g);
endfor
v=abs(v)^2
endfunction

valsByL = zeros(1,1000);

for k=1:1000,
valsByL(1,k)=sumSphCoeff(200*k,g);
endfor

plot( 1:1000, valsByL);```