;
; Uses the calibration residuals from Mike Mann's reconstructions to estimate
; the unresolved variance and its lag-1 autocorrelation.  These can then be
; used to estimate reconstruction errors at different time scales.  Does it
; for different fixed grids, to allow them to vary with time scale.
;
doar1=3    ; 0=no, 1=yes, 2=yes with fit to all lags, 3=1&2
doar2=3    ; 0=no, 1=yes, 2=yes with fit to all lags, 3=1&2
dolin=1    ; 0=no, 1=yes fit line to all lags (except 0)
;
fixedyr=[1000,1400,1600,1820]
nfix=n_elements(fixedyr)
;
nyr=1980-1902+1
rawdat=fltarr(4,nyr)
;
; Now prepare for plotting
;
loadct,39
multi_plot,nrow=4,ncol=2,layout='large'
if !d.name eq 'X' then begin
  wintim,ysize=750
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,/bold,font_size=12
endelse
;
nx=17
x=indgen(nx)
nxext=40
;
tsca=[1.,10,30,50,100,200]
nsca=n_elements(tsca)
;
allar1=fltarr(2,nfix)*!values.f_nan
allarn=fltarr(nx-1,2,nfix)*!values.f_nan
allarc=fltarr(2,2,nfix)*!values.f_nan
allarlin=fltarr(nxext,nfix)*!values.f_nan
allrlin=fltarr(nx-1,nfix)*!values.f_nan
;
allparam=fltarr(nxext+2,6,nfix)*!values.f_nan
;
allvif=fltarr(nsca,nfix)*!values.f_nan
allerr=fltarr(nsca,nfix)*!values.f_nan
;
for ifix = 0 , nfix-1 do begin
  ;
  ; Read residuals
  ;
  iyr=fixedyr[ifix]
  fn='nh-ad'+string(iyr,format='(I4)')+'-resid.dat'
  print
  print,'-------------------------------------'
  print,fn
  print,'-------------------------------------'
  openr,1,fn
  readf,1,rawdat
  close,1
  ;
  timey=reform(rawdat[0,*])
  resid=reform(rawdat[1,*])
  ;
  ; Plot residuals
  ;
  if ifix ne 0 then pause
  plot,timey,resid,title=fn,yrange=[-0.8,0.6],/ystyle
  oplot,!x.crange,[0,0],linestyle=1
  ;
  ; Analyse autoregressive and autocorrelation structure of residuals
  ;
  tjo_tsauto,resid,maxlag=nx-1,maxorder=7,wantorder=[1,2],$
    acf=acf,akaike=akaike,ar0=ar0,arwant=arwant,$
    bestorder=bestorder,/printakaike
  allparam[0:1,0,ifix]=ar0[*]        ; store AR(0) model
  allparam[0:2,1,ifix]=arwant[0,0:2] ; store AR(1) model
  allparam[0:3,2,ifix]=arwant[1,0:3] ; store AR(2) model
  allparam[0,3:5,ifix]=ar0[0]        ; store same mean in all other models
  ;
  xyouts,1905,!y.crange[0]+0.05,'SD='+string(ar0[1],format='(F5.3)')
  ;
  ; Plot autocorrelation function
  ;
  plot,x,acf,psym=def_sym(10),yrange=[-0.4,1],/ystyle,/xstyle,noclip=1
  oplot,!x.crange,[0,0],linestyle=1
  ;
  ; Compute and plot autocorrelation function for AR1 model
  ;
  resl1r=arwant[0,2]
  allar1[0,ifix]=resl1r
  acf_ar1=tjo_ar2acf(resl1r,maxlag=nx-1)
  if (doar1 mod 2) eq 1 then begin
    oplot,x,acf_ar1,thick=2
    xyouts,5,0.8,'o1='+string(resl1r,format='(F5.3)')
  endif
  ;
  ; Find best fit AR1 ACF to the actual ACF, then plot it
  ;
  iterfit,x,acf,0.,1.,resbest,rmse=rmse
  allar1[1,ifix]=resbest
  acf_ar1best=tjo_ar2acf(resbest,maxlag=nx-1)
  if doar1 ge 2 then begin
    oplot,x,acf_ar1best,thick=4
    xyouts,5,0.65,'o1fit='+string(resbest,format='(F5.3)')
  endif
  allparam[2,3,ifix]=resbest   ; store AR(1) best fit model
  ;
  ; For first two fixed grids, show the autocorrelation function for a
  ; second order model too
  ;
;  if ifix le 1 then begin
    ;
    ; Plot fitted AR2
    ;
    a=reform(arwant[1,2:3])
    acf_ar2=tjo_ar2acf(a,maxlag=nx-1)
    print
    print,'Standard AR2 parameters',a,format='(A,2F8.3)'
    if (doar2 mod 2) eq 1 then begin
      oplot,x,acf_ar2,thick=2,color=240
      xyouts,5,0.5,'o1,o2='+string(a,format='(2F6.3)')
    endif
    allarn[*,0,ifix]=acf_ar2[1:*]
    allarc[*,0,ifix]=a
    ;
    ; Iteratively fit AR2 to all of the autocorrelation function, choosing
    ; the best fit one
    ;
    if ifix le 1 then begin
    if doar2 ge 2 then begin
      ;
      ntry=10
      case ifix of
        0: begin
          amin=[0.23,0.23]  ; 0.00,0.00 : 0.07,0.00 : 0.17,0.10 : 0.19,0.19
          amax=[0.26,0.29]  ; 0.70,0.70 : 0.56,0.49 : 0.41,0.39 : 0.31,0.33
        end
        1: begin
          amin=[0.17,0.45]  ; 0.00,0.00 : 0.07,0.21 : 0.11,0.31 : 0.13,0.39
          amax=[0.21,0.49]  ; 0.70,0.70 : 0.42,0.56 : 0.32,0.56 : 0.26,0.53
        end
      endcase
      astep=(amax-amin)/float(ntry)
      aerr=fltarr(ntry+1,ntry+1)
      minerr=9999.
      mino1=9999.
      mino2=9999.
      miny=9999.
      for i = 0 , ntry do begin
        for j = 0 , ntry do begin
          a=[amin[0]+astep[0]*float(i),amin[1]+astep[1]*float(j)]
          acf_ar2best=tjo_ar2acf(a,maxlag=nx-1)
          aerr[i,j]=sqrt(total((acf-acf_ar2best)^2)/float(nyr))
          if aerr[i,j] lt minerr then begin
            minerr=aerr[i,j]
            mino1=a[0]
            mino2=a[1]
            miny=acf_ar2best
          endif
        endfor
      endfor
;      print,aerr,format='(11F7.4)'
;      print,minerr
      oplot,x,miny,thick=4,color=240
      xyouts,5,0.35,'o1,o2='+string([mino1,mino2],format='(2F6.3)')
      allarn[*,1,ifix]=miny[1:*]
      allarc[*,1,ifix]=[mino1,mino2]
      allparam[2:3,4,ifix]=[mino1,mino2]   ; store AR(2) best fit model
    endif
    endif
    ;
    ; Fit a line to the ACF
    ;
    xlin=x[1:nx-1]
    lincoeff=linfit(xlin,acf[1:nx-1])
    ylin=lincoeff[0]+lincoeff[1]*xlin
    allrlin[*,ifix]=ylin[*]
    if dolin ne 0 then oplot,xlin,ylin,color=120,thick=2
    ;
    ; Fit an AR(high-order) model to the linear fit to the ACF, extrapolated
    ; to zero and then extended with zeros
    ;
    ylinext=ylin
    if ylinext[nx-2] gt abs(lincoeff[1]) then begin
      ylinext=[ylinext,ylinext[nx-2]+lincoeff[1]]
    endif
    n2add=nxext-n_elements(ylinext)
    ylinext=[ylinext,replicate(0.,n2add)]
    ;
    arlin1=tjo_acf2ar(ylinext)
    allarlin[*,ifix]=arlin1[*]
    allparam[2:*,5,ifix]=arlin1[*]
    ;
    acf_arline=tjo_ar2acf(arlin1,maxlag=100)
;    plot,acf_arline,/xstyle
;    oplot,!x.crange,[0,0],linestyle=1
    ;
    print
    print,'Autoregressive parameters for straight-line fit'
    print,arlin1,format='(10F8.3)'
    ;
;  endif
  ;
  ; Compute reconstruction errors at different time scales, using
  ; both variance inflation factor formulae and by generating synthetic
  ; records.  For the latter, we take this opportunity to also estimate
  ; the noise term of each AR model.
  ;
  nrep=100L
  ngen=nyr*nrep
  allsyn=fltarr(ngen,3,2)
  print
  print,'Timescale-dependent error inflation factors, sqrt(v)'
  for i = 0 , nsca-1 do begin
    print,tsca[i],format='($,I4)'
    for j = 0 , 1 do begin
      r01=[1.,allar1[j,ifix]]
      ;
      v1=tjo_vif(r01)
      ;
      v2=tjo_vif(r01,avglen=tsca[i])
      if (ifix ge 2) and (j eq 1) then allvif[i,ifix]=v2
      ;
      ; Compare with syn
      if i eq 0 then begin
        synts=tjo_arsynth(r01[1],ngen)
        allsyn[*,0,j]=synts
        synts=reform(synts,nyr,nrep)
        varts=fltarr(nrep)
        for irep = 0 , nrep-1 do varts[irep]=variance(synts[*,irep])
        mvt=mean(varts)
        noisevar=ar0[1]^2 / mvt
        if j eq 0 then begin
          if abs(sqrt(noisevar)/allparam[1,1,ifix]-1.) gt 0.05 then $
            print,'CHECK NOISE!',sqrt(noisevar),allparam[1,1,ifix]
        endif else begin
          allparam[1,3,ifix]=sqrt(noisevar)
        endelse
      endif else begin
        synts=allsyn[*,0,j]
      endelse
      synts=reform(synts,ngen)
      ; Don't use "smooth" cos that only works for odd-length windows!
      nsmo=ngen-tsca[i]+1
      smots=fltarr(nsmo)
      for ismo = 0L , nsmo-1 do smots[ismo]=mean(synts[ismo:ismo+tsca[i]-1])
      v2syn=tsca[i]*variance(smots)/variance(synts)
      ;
      if finite(allarn[0,j,ifix]) eq 0 then begin
        v3=!values.f_nan
        v4=!values.f_nan
        v5=!values.f_nan
        v5syn=!values.f_nan
      endif else begin
        ;
        v3=tjo_vif(allarn[*,j,ifix])
        ;
        v4=tjo_vif(ar=allarc[*,j,ifix])
        ;
        v5=tjo_vif(ar=allarc[*,j,ifix],avglen=tsca[i])
        ;
        ; Compare with syn
        if i eq 0 then begin
          synts=tjo_arsynth(allarc[*,j,ifix],ngen)
          allsyn[*,1,j]=synts
          synts=reform(synts,nyr,nrep)
          varts=fltarr(nrep)
          for irep = 0 , nrep-1 do varts[irep]=variance(synts[*,irep])
          mvt=mean(varts)
          noisevar=ar0[1]^2 / mvt
          if j eq 0 then begin
            if abs(sqrt(noisevar)/allparam[1,2,ifix]-1.) gt 0.05 then $
              print,'CHECK NOISE!',sqrt(noisevar),allparam[1,2,ifix]
          endif else begin
            allparam[1,4,ifix]=sqrt(noisevar)
          endelse
        endif else begin
          synts=allsyn[*,1,j]
        endelse
        synts=reform(synts,ngen)
        ; Don't use "smooth" cos that only works for odd-length windows!
        nsmo=ngen-tsca[i]+1
        smots=fltarr(nsmo)
        for ismo = 0L , nsmo-1 do smots[ismo]=mean(synts[ismo:ismo+tsca[i]-1])
        v5syn=tsca[i]*variance(smots)/variance(synts)
        ;
      endelse
      ;
      print,sqrt([v1,v2,v2syn,v4,v5,v5syn]),format='($,6F7.3)'
      ;print,sqrt(v2),sqrt(v5),format='($,2F10.4)'
      ;print,sqrt(v3),sqrt(v4),format='($,2F10.4)'
;      print,sqrt(v2*(ar0[1]^2)/tsca[i]),sqrt(v5*(ar0[1]^2)/tsca[i]),format='($,2F10.4)'
    endfor
    ;
    ; Use the linear fit to the ACF
    v6=tjo_vif(allrlin[*,ifix],avglen=tsca[i])
    print,sqrt(v6),format='($,F7.3)'
    if ifix le 1 then allvif[i,ifix]=v6
    ;
    ; Compare with syn
    if i eq 0 then begin
      synts=tjo_arsynth(allarlin[*,ifix],ngen)
      allsyn[*,2,0]=synts
      synts=reform(synts,nyr,nrep)
      varts=fltarr(nrep)
      for irep = 0 , nrep-1 do varts[irep]=variance(synts[*,irep])
      mvt=mean(varts)
      noisevar=ar0[1]^2 / mvt
      allparam[1,5,ifix]=sqrt(noisevar)
    endif else begin
      synts=allsyn[*,2,0]
    endelse
    synts=reform(synts,ngen)
    ; Don't use "smooth" cos that only works for odd-length windows!
    nsmo=ngen-tsca[i]+1
    smots=fltarr(nsmo)
    for ismo = 0L , nsmo-1 do smots[ismo]=mean(synts[ismo:ismo+tsca[i]-1])
    v6syn=tsca[i]*variance(smots)/variance(synts)
    print,sqrt(v6syn),format='($,F7.3)'
    ;
    print
  endfor
  ;
  print
  print,' Timescale         VIF   sqrt(VIF)  standerr^2    standerr'
  for i = 0 , nsca-1 do begin
    vvv=allvif[i,ifix]
    eee=sqrt(vvv*allparam[1,0,ifix]^2/float(tsca[i]))
    allerr[i,ifix]=eee
    print,tsca[i],vvv,sqrt(vvv),eee^2,eee,format='(I10,4F12.4)'
  endfor
  ;
endfor
;
save,filename='mann_tderr.idlsave',$
  allvif,allerr,tsca,nsca,nfix,fixedyr
;
end