% Commands to perform generalized weighted least-squares regression.
% From the Wednesday, Sept 25th lecture in Math 456
% M.V.Wickerhauser


% Individual fitting functions.  Input v should be a column vector
f1 = @(v) ones (size (v));  % constant
f2 = @(v) exp (v-39);       % exponential, shifted 
f3 = @(v) cos (v);          % cosine

% Matrix whose three columns will be the fitting functions
F = @(v) [f1(v), f2(v), f3(v)];   % three column vectors

% Strike prices, equispaced from 38 to 41 by 0.50
ks = (38:0.50:41)'  % transpose to get a column vector

% Call premiums by strike price:
cp = [1.25 0.90 0.44 0.21 0.08 0.03 0.01]'  % transpose to column vector

% Open interest vector injects into the diagonal Weight matrix:
oi = [710 609 989 5360 10063 5717 7770];  % OK to have a row vector here
W = diag(oi)

% Find the expansion coefficients in f1,f2,f3 by generalized weighted
%   least-squares regression:
pW = (F(ks)'*W*F(ks))\(F(ks)'*W*cp)


% Dense x points for plotting,
x = 38:0.02:41;  x=x'; % transpose to get a column vector

% Spline interpolation of the original data:
sfit = spline(ks,cp); ysi = ppval(sfit, x);

% Weighted regression function:
ywr = F(x)*pW;

% Plot the points, the spline interpolation, and the weighted
%  regression all on one graph:
plot(ks,cp,"b+", x,ysi, "b-", x,ywr, "k--");