plmt2resid.m

function varargout=plmt2resid(plmt,thedates,fitwhat,givenerrors)
% [ESTresid,thedates,ESTsignal,rmset,rmsst,varet,varst]
%             =PLMT2RESID(plmt,thedates,fitwhat,givenerrors)
%
% Takes a time series of spherical harmonic coefficients and fits a desired
% combination of functions (e.g. secular, annual, semiannual, etc.)
%
% INPUT:
%
% plmt        The time series of spherical harmonic coefficients.  This
%               should be a three dimensional matrix (not a cell array), 
%               where the first dimension is time, and dimensions 2 and 3 
%               are as in the traditional lmcosi matrix.
% thedates    An array of dates corresponding to the plm thedateseries.  These
%               should be in Matlab's date format. (see DATENUM)
% fitwhat     The functions that you would like to fit to the time series
%               data.  The format of this array is as follows:
%               [mean secular periodic1 periodic2 periodic3 etc...] where
%               - mean is either 0/1 to turn off/on this function
%               - secular is either 0/1 to turn off/on this function
%               - periodic1 is the period in days of a function (i.e. 365.0)
%              Any # of desired periodic functions can be included. 
% givenerrors  These are given errors, if you have them.  In this case a
%                weighted inversion is performed.  givenerrors should be the
%                same dimensions of plmt.
%            
% OUTPUT:
%
% ESTresid    Residual time series for each pair of cos/sin coefficients
%             [nmonths x Ldata x 4] with a hyper-lmcosi format
%             ESTresid(:,:,1) the degree
%             ESTresid(:,:,2) the order
%             ESTresid(:,:,3) the cosine coefficient
%             ESTresid(:,:,4) the sine coefficient
% thedates    The dates that belong to these time series
% ESTsignal   The least-squares fitted function for each cos/sin
%             coefficient evaluated at those months in the same format 
% rmset       The root mean square of the estimated residual coefficients
% rmsst       The root mean square of the estimated signal coefficients
% varet       The variance of the estimated residual coefficients
% varst       The variance of the estimated signal coefficients
%
% SEE ALSO:
%
% Last modified by charig-at-princeton.edu 5/29/2023
% Last modified by fjsimons-at-alum.mit.edu 5/26/2011

defval('xver',0)
defval('plmt',fullfile(getenv('IFILES'),'GRACE','CSR_RL05_alldata.mat'))
warning('off','MATLAB:divideByZero');

if isstr(plmt)
  % Load the file specified
  a=load(plmt);
  % These are the potential coefficients
  plmt=a.potcoffs(:,:,1:4);
  % These would be with the "formal" errors
  % givenerrros=a.potcoffs(:,:,[1 2 5:6]);
  % Much better to use the "calibrated" errors
  %givenerrors=a.cal_errors(:,:,1:4);
  % Now get the thedates for each of those months
  thedates=a.thedates;
  % Get rid of the rest
  clear a
end

% Initialize/Preallocate
defval('givenerrors',ones(size(plmt)));
defval('fitwhat',[3 181.0 365.0])
defval('P2ftest',0);
defval('P3ftest',0);

% Round to 1 January of the first available and next last available year 
bounds=[datevec(thedates(1)); datevec(thedates(end))];
%x0 = datenum(bounds(1)  ,1,1);
%x1 = datenum(bounds(2)+1,1,1);


[i,j,k] = size(plmt);

% Initialize the residuals
ESTresid =zeros(size(plmt));

% Initialize the evaluated fitted function set
ESTsignal=zeros(size(plmt));

% Put in the degrees and the orders
ESTresid(:,:,1:2) =plmt(:,:,1:2);
ESTsignal(:,:,1:2)=plmt(:,:,1:2);

% Figure out the orders and degrees of this setup
Ldata=addmup(j,'r');
[dems,dels]=addmon(Ldata);
difer(squeeze(plmt(1,:,2))-dems(:)',[],[],NaN)
difer(squeeze(plmt(1,:,1))-dels(:)',[],[],NaN)

% How many data?
nmonths=length(thedates);

% We will do a scaling to improve the solution
mu1 = mean(thedates); % mean
mu2 = std(thedates); % standard deviation
% Make a new x-vector with this information
xprime = (thedates - mu1)/mu2;

% The frequencies being fitted in [1/days]
omega = 1./[fitwhat(2:end)];
% Rescale these to our new xprime
omega = omega*mu2;

% How many periodic components?
lomega=length(omega);


%%%
% G matrix assembly
%%%

% We will have the same number of G matrices as order of polynomial fit.
% These matrices are smallish, so make all 3 regardless of whether you want
% them all.
G1 = []; % For line fits
G2 = []; % For quadratic fits
G3 = []; % For cubic fits
% Mean term
if fitwhat(1) >= 0
  G1 = [G1 ones(size(xprime'))];
  G2 = [G2 ones(size(xprime'))];
  G3 = [G3 ones(size(xprime'))];
end
% Secular term
if fitwhat(1) >= 1
  G1 = [G1 (xprime)'];
  G2 = [G2 (xprime)'];
  G3 = [G3 (xprime)'];
end
% Quadratic term
if fitwhat(1) >= 2
  G2 = [G2 (xprime)'.^2];
  G3 = [G3 (xprime)'.^2];
end
% Cubic term
if fitwhat(1) == 3
  G3 = [G3 (xprime)'.^3];
end
% Periodic terms
if ~isempty(omega)
  % Angular frequency in radians/(rescaled day) of the periodic terms
  th_o= repmat(omega,nmonths,1)*2*pi.*repmat((xprime)',1,lomega);
  G1 = [G1 cos(th_o) sin(th_o)];
  G2 = [G2 cos(th_o) sin(th_o)];
  G3 = [G3 cos(th_o) sin(th_o)];

end % end periodic

% Initilization Complete


%%%
% Solving
%%%

% Loop over coefficients starting from the (1,0) terms and do fitting
for index=2:j
    % If we have a priori error information, create a weighting matrix, and
    % change the G and d matrices to reflect this.  Since each coefficient
    % has its own weighting, we have to invert them separately.
    W3 = diag([1./squeeze(givenerrors(:,index,3))]);
    % Isolate the data over time
    d = squeeze(plmt(:,index,3:4));
    G1wc = W3*G1;
    G2wc = W3*G2;
    G3wc = W3*G3;

    if dems(index)==0
      % The sine coefficient and errors are 0. So prevent an error.
      G1ws =  G1;
      G2ws =  G2;
      G3ws =  G3;
      dw = [W3*d(:,1) d(:,2)];
    else
      W4 = diag([1./squeeze(givenerrors(:,index,4))]);
      G1ws =  W4*G1;
      G2ws =  W4*G2;
      G3ws =  W4*G3;
      dw = [W3*d(:,1) W4*d(:,2)];
    end

    % Special case from slept2resid that we don't use here
    th=th_o;
    
    %%%
    % First order polynomial
    %%%

    % Do the fitting by minimizing least squares which is exactly right
    mL2_1 = [(G1wc'*G1wc)\(G1wc'*dw(:,1)) (G1ws'*G1ws)\(G1ws'*dw(:,2))];
    % mL2 is m by 2, where col 1 is for the C terms, col 2 for S terms
    
    % Use the model parameters to make the periodic amplitude in time
    startP = length(mL2_1) - 2*lomega + 1;
    Camp = [mL2_1(startP:(startP+lomega-1),1) mL2_1((startP+lomega):end,1)];
    Camp = sqrt(Camp(:,1).^2 + Camp(:,2).^2);
    
    Samp = [mL2_1(startP:(startP+lomega-1),2) mL2_1((startP+lomega):end,2)];
    Samp = sqrt(Samp(:,1).^2 + Samp(:,2).^2);
    
    % Use the model parameters to make the periodic phase in time
    Cphase = [mL2_1(startP:(startP+lomega-1),1) mL2_1((startP+lomega):end,1)];
    Cphase = atan2(Cphase(:,1),Cphase(:,2));

    Sphase = [mL2_1(startP:(startP+lomega-1),2) mL2_1((startP+lomega):end,2)];
    Sphase = atan2(Sphase(:,1),Sphase(:,2));
    
    % Assemble the estimated signal function, evaluated at 'thedates'
    fitfnC1 = 0;
    fitfnS1 = 0;
    
    % Start adding things
    fitfnC1 = fitfnC1 + mL2_1(1,1) + mL2_1(2,1)*(xprime);
    fitfnS1 = fitfnS1 + mL2_1(1,2) + mL2_1(2,2)*(xprime);
    
    if ~isempty(omega)
    % Add the sum over all the periodic components periodics
    fitfnC1 = fitfnC1 + ...
             sum(repmat(Camp,1,nmonths).*sin(th'+repmat(Cphase,1,nmonths)),1);
    fitfnS1 = fitfnS1 + ...
             sum(repmat(Samp,1,nmonths).*sin(th'+repmat(Sphase,1,nmonths)),1);
    end
    
    % This is completely equivalent with the above for the purpose of residual
    % determination
    if xver==1
      % sum(repmat(mL2(startP:(startP+lomega-1),1),1,nmonths).*cos(th'),1)+...
      %       sum(repmat(mL2((startP+lomega):end,1),1,nmonths).*sin(th'),1)
      % sum(repmat(mL2(startP:(startP+lomega-1),2),1,nmonths).*cos(th'),1)+...
      %       sum(repmat(mL2((startP+lomega):end,2),1,nmonths).*sin(th'),1)
      clf
      subplot(211)
      plot(thedates,d(:,1),'+')
      hold on
      plot(thedates,fitfnC,'r');
      axis tight
      datetick('x',28)
      title(sprintf('cosine component at l = %i and m = %i',...
                    dels(index),dems(index)))
      
      subplot(212)
      plot(thedates,d(:,2),'+')
      hold
      plot(thedates,fitfnS,'r');
      axis tight
      datetick('x',28)
      title(sprintf('sine component at l = %i and m = %i',...
                    dels(index),dems(index)))
      pause(2)
    end
    
    % Compute the residual time series for this pair of cosine and sine
    % coefficients 
    resid1 = d - [fitfnC1' fitfnS1'];
    
    % Here's the definition of the residual at lm vs time
    ESTresid(:,index,3:4) = resid1;
    % Here's the definition of the signal at lm vs time
    ESTsignal(:,index,3:4) = [fitfnC1' fitfnS1'];

    % Get the residual sum of squares for later F tests
    rss1 = sum(resid1.^2,1);



    %%%
    % Second order polynomial
    %%%

    % Now repeat that all again with second order polynomial, if we want
    if fitwhat(1) >= 2  
        % Do the fitting by minimizing least squares which is exactly right
        mL2_2 = [(G2wc'*G2wc)\(G2wc'*dw(:,1)) (G2ws'*G2ws)\(G2ws'*dw(:,2))];
      

        % Use the model parameters to make the periodic amplitude in time
        startP = length(mL2_2) - 2*lomega + 1;
        Camp = [mL2_2(startP:(startP+lomega-1),1) mL2_2((startP+lomega):end,1)];
        Camp = sqrt(Camp(:,1).^2 + Camp(:,2).^2);
    
      Samp = [mL2_2(startP:(startP+lomega-1),2) mL2_2((startP+lomega):end,2)];
      Samp = sqrt(Samp(:,1).^2 + Samp(:,2).^2);
    
      % Use the model parameters to make the periodic phase in time
      Cphase = [mL2_2(startP:(startP+lomega-1),1) mL2_2((startP+lomega):end,1)];
      Cphase = atan2(Cphase(:,1),Cphase(:,2));

      Sphase = [mL2_2(startP:(startP+lomega-1),2) mL2_2((startP+lomega):end,2)];
      Sphase = atan2(Sphase(:,1),Sphase(:,2));


      % Assemble the estimated signal function, evaluated at 'thedates'
      fitfnC2 = 0;
      fitfnS2 = 0;
    
      % Start adding things
      fitfnC2 = fitfnC2 + mL2_2(1,1) + mL2_2(2,1)*(xprime) + mL2_2(3,1)*(xprime).^2;
      fitfnS2 = fitfnS2 + mL2_2(1,2) + mL2_2(2,2)*(xprime) + mL2_2(3,2)*(xprime).^2;
    
      if ~isempty(omega)
         % Add the sum over all the periodic components periodics
         fitfnC2 = fitfnC2 + ...
             sum(repmat(Camp,1,nmonths).*sin(th'+repmat(Cphase,1,nmonths)),1);
         fitfnS2 = fitfnS2 + ...
             sum(repmat(Samp,1,nmonths).*sin(th'+repmat(Sphase,1,nmonths)),1);
      end

    
      % Compute the residual time series for this pair of cosine and sine
      % coefficients 
      resid2 = d - [fitfnC2' fitfnS2'];
   
      
      % Get the residual sum of squares
      rss2 = sum(resid2.^2,1);
      % Calculate an F-score for this new fit
      fratioP2C = (rss1(1) - rss2(1))/1/(rss2(1)/(nmonths-length(mL2_2(:,1))));
      fratioP2S = (rss1(2) - rss2(2))/1/(rss2(2)/(nmonths-length(mL2_2(:,2))));
      fscore2C = finv(.95,1,nmonths-length(mL2_2(:,1)));
      fscore2S = finv(.95,1,nmonths-length(mL2_2(:,2)));
      if fratioP2C > fscore2C
         P2Cftest = 1;
         % We pass, so update the signal and residuals with this new fit
         ESTresid(:,index,3) = resid2(:,1);
         ESTsignal(:,index,3) = fitfnC2';
      else
         P2ftest = 0;
      end
      if fratioP2S > fscore2S
         P2ftest = 1;
         % We pass, so update the signal and residuals with this new fit
         ESTresid(:,index,4) = resid2(:,2);
         ESTsignal(:,index,4) = fitfnS2';
      else
         P2ftest = 0;
      end

    end % end second order
    


    %%%
    % Third order polynomial
    %%%
    
    % Now repeat that all again with third order polynomial, if we want
    if fitwhat(1) >= 3
        % Do the fitting by minimizing least squares which is exactly right
        mL2_3 = [(G3wc'*G3wc)\(G3wc'*dw(:,1)) (G3ws'*G3ws)\(G3ws'*dw(:,2))];
      

        % Use the model parameters to make the periodic amplitude in time
        startP = length(mL2_3) - 2*lomega + 1;
        Camp = [mL2_3(startP:(startP+lomega-1),1) mL2_3((startP+lomega):end,1)];
        Camp = sqrt(Camp(:,1).^2 + Camp(:,2).^2);
    
      Samp = [mL2_3(startP:(startP+lomega-1),2) mL2_3((startP+lomega):end,2)];
      Samp = sqrt(Samp(:,1).^2 + Samp(:,2).^2);
    
      % Use the model parameters to make the periodic phase in time
      Cphase = [mL2_3(startP:(startP+lomega-1),1) mL2_3((startP+lomega):end,1)];
      Cphase = atan2(Cphase(:,1),Cphase(:,2));

      Sphase = [mL2_3(startP:(startP+lomega-1),2) mL2_3((startP+lomega):end,2)];
      Sphase = atan2(Sphase(:,1),Sphase(:,2));


      % Assemble the estimated signal function, evaluated at 'thedates'
      fitfnC3 = 0;
      fitfnS3 = 0;
    
      % Start adding things
      fitfnC3 = fitfnC3 + mL2_3(1,1) + mL2_3(2,1)*(xprime) + mL2_3(3,1)*(xprime).^2;
      fitfnS3 = fitfnS3 + mL2_3(1,2) + mL2_3(2,2)*(xprime) + mL2_3(3,2)*(xprime).^2;
    
      if ~isempty(omega)
         % Add the sum over all the periodic components periodics
         fitfnC3 = fitfnC3 + ...
             sum(repmat(Camp,1,nmonths).*sin(th'+repmat(Cphase,1,nmonths)),1);
         fitfnS3 = fitfnS3 + ...
             sum(repmat(Samp,1,nmonths).*sin(th'+repmat(Sphase,1,nmonths)),1);
      end

    
      % Compute the residual time series for this pair of cosine and sine
      % coefficients 
      resid3 = d - [fitfnC3' fitfnS3'];


      % Get the residual sum of squares
      rss3 = sum(resid3.^2,1);
      % Calculate an F-score for this new fit
      fratioP3C = (rss1(1) - rss3(1))/1/(rss3(1)/(nmonths-length(mL2_3(:,1))));
      fratioP3S = (rss1(2) - rss3(2))/1/(rss3(2)/(nmonths-length(mL2_3(:,2))));
      fscore3C = finv(.95,1,nmonths-length(mL2_3(:,1)));
      fscore3S = finv(.95,1,nmonths-length(mL2_3(:,2)));
      if fratioP3C > fscore3C
         P3Cftest = 1;
         % We pass, so update the signal and residuals with this new fit
         ESTresid(:,index,3) = resid3(:,1);
         ESTsignal(:,index,3) = fitfnC3';
      else
         P3ftest = 0;
      end
      if fratioP3S > fscore3S
         P3ftest = 1;
         % We pass, so update the signal and residuals with this new fit
         ESTresid(:,index,4) = resid3(:,2);
         ESTsignal(:,index,4) = fitfnS3';
      else
         P3ftest = 0;
      end

      
    end % end third order
    
    % Make the matrix ftests
    ftests(index,:) = [0 P2ftest P3ftest];


end


% Compute the RMS of the estimated signal
rmsst = squeeze(ESTsignal(1,:,1:2));
rmsst(:,3:4) = squeeze( sqrt( mean( ESTsignal(:,:,3:4).^2,1 ) ) );

% Compute the RMS of the estimated residual
rmset = squeeze(ESTresid(1,:,1:2));
rmset(:,3:4) = squeeze( sqrt( mean( ESTresid(:,:,3:4).^2,1 ) ) );

% Compute the variance of the estimated signal
varst = squeeze(ESTsignal(1,:,1:2));
varst(:,3:4) = squeeze( var( ESTsignal(:,:,3:4),1 ) );

% Compute the variance of the estimated signal
varet = squeeze(ESTresid(1,:,1:2));
varet(:,3:4) = squeeze( var( ESTresid(:,:,3:4),1 ) );

% Collect output
varns={ESTresid,thedates,ESTsignal,rmset,rmsst,varet,varst};
varargout=varns(1:nargout);