;
; Computes regressions on full, high and low pass MEAN timeseries of MXD
; anomalies against full NH temperatures.
; THIS IS FOR THE AGE-BANDED (ALL BANDS) STUFF OF HARRY'S
;
; Specify period over which to compute the regressions (stop in 1940 to avoid
; the decline
;
perst=1881.
peren=1960.
;
openw,1,'recon1_nhemi.out'
thalf=10.
;
; Read in hemispheric timseries and filter
;
restore,'instrts_nh.idlsave'
totval=total(nhts(*,3:8),2,/nan)
totn=total(finite(nhts(*,3:8)),2)
keeplist=where(totn gt 2,nkeep)
lvx=x
lvy=totval*!values.f_nan
lvy(keeplist)=totval(keeplist)/float(totn(keeplist)) 
filter_cru,thalf,tsin=lvy,tslow=lvylow,tshigh=lvyhi,/nan
;
; Read in MXD timeseries and filter
;
restore,filename='densadj_all(330).idlsave'
if total(abs(lvx-x)) ne 0 then message,'Incompatible years'
y1=densall
filter_cru,thalf,tsin=y1,tslow=ylow1,tshigh=yhi1,/nan
;
; Now correlate and regress them
;
printf,1,'Correlations and regression coefficients between normalised timeseries'
keeplist=where(finite(y1+lvy) and (x ge perst) and (x le peren),nkeep)
x=x(keeplist)
ts1=y1(keeplist)
ts2=lvy(keeplist)
r1=correlate(ts1,ts2)
dum=linfit(ts1,ts2,sigma=err_sigma)
c1=dum
printf,1,'Full timeseries:',r1,c1
;
tsa=ts1(0:nkeep-2)
tsb=ts1(1:nkeep-1)
;mxd_ar1=correlate(tsa,tsb)
mxd_ar1=a_correlate(ts1,[-1])
tsa=ts2(0:nkeep-2)
tsb=ts2(1:nkeep-1)
;nht_ar1=correlate(tsa,tsb)
nht_ar1=a_correlate(ts2,[-1])
printf,1,'AR1 for MXD and NHEMI:',mxd_ar1,nht_ar1
;
multi_plot,nrow=2
if !d.name eq 'X' then window,ysize=800
tspred=c1(0)+c1(1)*ts1
tswant=ts2
plot,tspred,ts2,psym=1,$
  xtitle='Scaled age-banded MXD anomaly  (!Uo!NC)',$
  ytitle='Northern Hemisphere temperature anomaly  (!Uo!NC)',$
  /xstyle,xrange=[-0.6,0.3],$
  /ystyle,yrange=[-0.6,0.3]
oplot,!x.crange,[0.,0.],linestyle=1
oplot,[0.,0.],!y.crange,linestyle=1
dum=linfit(tspred,ts2,sigma=se)
oplot,!x.crange,!x.crange*dum(1)+dum(0)
oplot,!x.crange,!x.crange*dum(1)+dum(0)+se(0)
oplot,!x.crange,!x.crange*dum(1)+dum(0)-se(0)
oplot,!x.crange,!x.crange*(dum(1)+se(1))+dum(0)
oplot,!x.crange,!x.crange*(dum(1)-se(1))+dum(0)
oplot,!x.crange,!x.crange*(dum(1)+se(1))+dum(0)+se(0)
oplot,!x.crange,!x.crange*(dum(1)-se(1))+dum(0)-se(0)
oplot,!x.crange,!x.crange*(dum(1)+se(1))+dum(0)-se(0)
oplot,!x.crange,!x.crange*(dum(1)-se(1))+dum(0)+se(0)
;
ts1=yhi1(keeplist)
ts2=lvyhi(keeplist)
r2=correlate(ts1,ts2)
dum=linfit(ts1,ts2)
c2=dum
printf,1,'High-pass      :',r2,c2
;
ts1=ylow1(keeplist)
ts2=lvylow(keeplist)
r3=correlate(ts1,ts2)
dum=linfit(ts1,ts2)
c3=dum
printf,1,'Low-pass       :',r3,c3
;
; Now compute the rms error between the reconstruction and the original
; timeseries
;
tserr=tswant-tspred
rmserr=sqrt( total( tserr^2 ) / float(n_elements(tserr)) )
printf,1,'RMS error between NH temperature and age-banded-MXD reconstruction'
printf,1,rmserr
printf,1,'Uncertainties surrounding regression coefficients'
printf,1,err_sigma
;
; Now do some example reconstructions (also for 1912 & 1884 we know the
; instrumental values)
;
printf,1,' '
ele1=where(x eq 1884)
printf,1,'Instrumental NH temp for 1884 is',tswant(ele1(0))
ele1=where(x eq 1912)
printf,1,'Instrumental NH temp for 1912 is',tswant(ele1(0))
printf,1,' '
;
exyr=[1601,1816,1641,1453,1912,1884]
exval=[-6.9026,-4.3256,-4.3124,-4.2402,-3.3282,-2.8874]
nex=n_elements(exyr)
printf,1,'Reconstructions with error bounds'
for i = 0 , nex-1 do begin
  xx=exval(i)
  y1=xx*(c1(1)+err_sigma(1))+c1(0)-err_sigma(0)-rmserr
  y2=xx*c1(1)+c1(0)
  y3=xx*(c1(1)-err_sigma(1))+c1(0)+err_sigma(0)+rmserr
  printf,1,exyr(i),y1,y2,y3
endfor
;
printf,1,' '
printf,1,'Computations carried out over the period ',perst,peren
printf,1,' '
printf,1,'To separate low and high frequency components, a gaussian weighted'
printf,1,'filter was used with a half-width (years) of ',thalf
;
close,1
;
end