pro mkcalibratemulti,por,pand,fitcoeff=fitcoeff,autocoeff=autocoeff,$
  nlag=nlag
;
; Given a set of predictor time series (por[n,t]) and a predictand time series
; (pand[t]), with no missing data, multiple linear (least-squares)
; regression is performed to calibrate the predictors.
; fitcoeff is returned containing {r,alpha,beta,SEalpha,SEbeta,sig,rms},
; where r=correlation(prediction,pand), alpha and beta are the offset and
; scaling coefficients to be applied to the predictors, their standard errors
; are also given, the significance (as a % confidence) is returned (taking
; into account autocorrelation), and also the root-mean-squared error after
; calibration.  autocoeff is returned containing {lag1r of por, lag1r of pand,
; and lag1r of residuals, effective independent sample size}.
;
; When computing the signficance of the correlation, the number of lagged
; autocorrelations to use can be specified in nlag (default is 1).  This
; is usually only needed when using smoothed series that no longer look
; like simple red noise (when nlag=9 typically works well).
; The lag1r of each predictor is computed, and then only the maximum of these
; is used (to reduce sample size), although all are returned.
;
; NB. The multiple regression function does not compute a standard error
; for the constant (intercept) coefficient, so this is set to missing.
;
; Check time series lengths and number of predictors
;
psize=size(por)
case psize(0) of
  1: begin
    nval=psize(1)
    npred=1
    por=reform(por,npred,nval)
  end
  2: begin
    nval=psize(2)
    npred=psize(1)
  end
  else: message,'por[n,t] must be 1D or 2D'
endcase
if n_elements(pand) ne nval then message,'Incompatible time series'
if n_elements(nlag) eq 0 then nlag=1
;print,npred,nval
;
; Compute multiple regression
;
wts=fltarr(nval)+1.
bcoeff=regress(por,pand,wts,/relative_weight,predval,acoeff,bse,fvalue,$
  sepr,mulr)
bcoeff=reform(bcoeff)
;print,bcoeff
;print,sepr
;print,fvalue
;
; Compute residuals
resids=pand-predval
;
; Compute RMS error
sqerr=resids^2
rmse=sqrt( total(sqerr)/nval )
;
; Compute lag-1 autocorrelations
alag=indgen(nlag)+1
ac1all=fltarr(npred,nlag)
for ipred = 0 , npred-1 do begin
  ac1all(ipred,*)=a_correlate(reform(por(ipred,*)),alag)
endfor
dummy=max(reform(ac1all(*,0)),maxele)
;print,dummy
ac1por=reform(ac1all(maxele,*))
ac1pand=a_correlate(pand,alag)
ac1resid=a_correlate(resids,alag)
;
; Compute significance of r(por,pand), taking autocorrelation into account
rr = (ac1por(0) > 0.) * (ac1pand(0) > 0.)
if nlag eq 1 then ess=(nval-2.)*(1.-rr)/(1.+rr)+2. $
             else ess=(nval-2.)/(1.+2.*total(rr))+2.
v1=npred
v2=(ess-npred-1) > 1
tsig=f_pdf(fvalue,v1,v2)
tsig=tsig*100.
;
fitcoeff = { r: mulr, alpha: acoeff, beta: bcoeff, sealpha: !values.f_nan,$
  sebeta: bse, sig: tsig, rms: rmse }
autocoeff= { ac1por: reform(ac1all(*,0)), ac1pand: ac1pand(0), $
  ac1resid: ac1resid(0), ess:ess }
;
return
end