function [muR,SEmuR,PvR,SvR] = remltrend( data )
% [muR,SEmuR,PvR,SvR] = remltrend( data )
%  function to estimate trend and Process variance using methods
%  presented in Staples et al
%
%  NOTE:  this function calls optimize.m function by Staples et al.
%  
%  muR, SEmuR, PvR, SvR are trend, standard error of trend,
%  process variance, non-process variance estimates
%  from Staples et al method
%
%  if either PvR or SvR are equal to 1.0000e-005, these are boundary
%    estimates for the parameters;  see Staples et al
%
%  data is a column vector of population estimates!!

    [years,rdata]=size(data);

%  w is vector of log differeces in population size
    w = log(data(2:years)) - log(data(1:years-1));
    [wr wc]= size(w);
    
    % construct initial estimates of V - cov matrix for w
    %  using method of moments
    sig1hat = var(w);
    covw = cov(w(2:wr),w(1:wr-1));
    sig2hat = covw(1,2);
    sigs = [sig1hat;sig2hat];
    V_hat = diag(sig1hat*ones(wr,1)) + diag(sig2hat*ones(wr-1,1),1)...
        + diag(sig2hat*ones(wr-1,1),-1);
    
    
    %  generate U matrix such that U'*X = 0 
    %  This annhilates fixed effect (mean) for REML estimation
    U = diag(-ones(wr,1)) + diag(ones(wr-1,1),-1);
    U = U(:,1:wr-1);
    
    % create time-series (z) with fixed mean annhilated
    z = U'*w;
    [zr zc]=size(z);
    
    %  conditions for optimization routing -- see MATLAB fmincon function
    %  necessary to restrict the parameter space such that V is invertable
    A = eye(2);
    A(1,1) = -1;
    A(3,:)=[-1 -2];
    A(4,:)=[-1 2];
    b = [0;0;-.00001;-.00001];
    
    %  REML procedure
    options = optimset('LargeScale','off','Display','off',...
        'MaxFunEvals',200,'MaxIter',200);
    %fmincon(FUN,XO,A,B,Aeq,Beq,LB,UB,NONLCON, OPTIONS, P1, P2...)
    x = fmincon('optimize',sigs,A,b,[],[],[0;-2*sig1hat],[2*sig1hat;0],[],options,z,U);
    
    %  model matrix
    X = ones(wr,1);
        
    %  REML estimates for variance paramters
    sig1 = x(1);
    sig2 = x(2);
    
    %  estimate of process variance
    k=[1 2];
    PvR = k*x;
    
    %  estimate of sampling variance
    SvR = -sig2;  
    
    %  calculate Variance for trend estimate
    c2 = sig2*ones(wr-1,1);
    V = diag(sig1*ones(wr,1)) + diag(c2,1) + diag(c2,-1);
    iV = inv(V);
    phi = inv(X'*iV*X);  % var(trend estimate)
   
    %   estimate of trend and SE trend
    muR = phi*X'*iV*w ;  
    SEmuR = sqrt(phi);