Added iWLS scripts from Homer2

function [dmoco, beta, tstat, pval, sigma, CovB, dfe, w, P, f] = ar_glm_final( d,X,tune )

    % preallocation

    dmoco = zeros(size(d));

    beta = zeros([size(X,2) size(d,2)]);

    tstat = beta;

    pval = beta;

    w = zeros(size(d));

    Sall=[];

    for j = 1:size(d,2)

        y = d(:,j);

        B = X \ y;

        B0 = 1e6*ones(size(B));

        iter = 0;

        maxiter = 10;

        while norm(B-B0)/norm(B0) > 1e-2 && iter < maxiter 

            B0 = B;

            res = y-X*B;

            a = robust_ar_fit(res);

            f = [1; -a(2:end)];

            Xf = filter(f,1,X);

            yf = filter(f,1,y);

            [B, S] = robustfit(Xf,yf,[],[],'off');

            iter = iter + 1;

        end

        % moco data & statistics

        w(:,j) = S.w;

%         dmoco(:,j) = ymoco;

        beta(:,j) = B;

        tstat(:,j) = B./sqrt(diag(S.covb));

        pval(:,j) = 2*tcdf(-abs(tstat(:,j)),S.dfe);

        P(j,1) = length(a)-1;

        sigma(:,j)=S.s;

        CovB(:,:,j)=S.covb;

        dfe=S.dfe;

    end

   dmoco=filter(f,1,d); 

end

% function w = biweight( y,c )

%     sig = mad(y,1)/.6745;

%     y = y/sig;

%     

%     lst = y < c & y > -c;

%     w = zeros( size(y) );

%     w(lst) = abs(1 - (y(lst)/c).^2);

% end

% function F = getFilterMatrix( a,ly )

%     P = length(a)-1;

%     F = speye(ly);

%     F = spdiags(repmat(-a(2:end)',[ly 1]),-1:-1:-P,F);

% end

function [beta, P] = robust_ar_fit( y )

	warning('off','stats:statrobustfit:IterationLimit')

    yf = y;

    yb = flipud(y);

    bic0 = 1; bic = 0;

    P = 0; Pmax = 100;

    beta0 = []; beta = [];

    while bic0 > bic && P < Pmax

        beta0 = beta;

        bic0 = bic;

        P = P+1; 

        % forward design matrix

        Xf = zeros(length(yf)-P,P+1);

        Xf(:,1) = ones(length(yf)-P,1);

        for j = 1:P

            Xf(:,j+1) = yf(j+1:end-P+j);

        end

        % backward design matrix

        Xb = zeros(length(yb)-P,P+1);

        Xb(:,1) = ones(length(yb)-P,1);

        for j = 1:P

            Xb(:,j+1) = yb(j+1:end-P+j);

        end

        %  weighted least squares

        y = [yf(1:end-P); yb(1:end-P)];

        X = [Xf; Xb];

        beta = pinv(X) * y;

        r = y - X*beta;

        w = biweight(r);

        beta = pinv((repmat(w,[1 size(X,2)]).*X))* (w.*y);

        r = y - X*beta;

%         % robust regression

%         y = [yf(1:end-P); yb(1:end-P)];

%         X = [Xf; Xb];

%         [beta, ~] = robustfit(X,y,[],[],'off');

%         r = y-X*beta;

        % AIC & BIC

        n = length(r);

        LL = -n/2*log(2*pi*var(r))  - n/2;

        bic = -LL + P/2*log(n);

    end

    if P > 1

        beta = beta0;

    end

end

function w = biweight( r )

    c = 4.685;

    sig = mad(r,1)/.6745;

    r = r/sig;

    lst = r < c & r > -c;

    w = zeros( size(r) );

    w(lst) = abs(1 - (r(lst)/c).^2);

end