% For Part 4b by Stig Mortensen
% use linear dependence between y(t) and y(t-1)
% least squares solution -  common mu

p =  % number of points to evaluate on each axis
xp = % vector of p points for x axis
y  = % vector of p points for y axis 
h  = % fixed bandwidth
order = 2;
X = []; W=[]; Yt = [];
for i = 1:length(xp)
    idx = find(data(:,1)<xp(i)+h/2 & data(:,1)>=xp(i)-h/2);
    n = length(idx);
    Yt = [Yt ; data(idx,3)];
    Yt_1 = data(idx,2);
    Xt_1 = data(idx,1)-xp(i);
    switch order
        case 2
            X = [ X ; ...
                ones(n,1) zeros(n,(i-1)*3) Yt_1  Yt_1.*Xt_1  Yt_1.*(Xt_1.^2) ...
                zeros(n,3*(p-i))];
        case 1
            X = [ X ; ...
                ones(n,1) zeros(n,(i-1)*2) Yt_1  Yt_1.*Xt_1  ...
                zeros(n,2*(p-i))];
    end
    u = abs(Xt_1)/max(abs(Xt_1)); %distance to xp scaled to [0 1]
    W = [W ;(1-u.^3).^3; ];
end
%W = diag(W);
W = spdiags(W,0,length(W),length(W));
tic,theta = (X'*W*X)\(X'*W*Yt); ls = toc
for i=1:p
    ZZ(:,i) = theta(1) + theta(2+(i-1)*(order+1))*y;
end