# -*- coding: utf-8 -*-
"""
Created on Mon Nov  8 14:05:29 2021

https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.linregress.html
https://hpklima.blogspot.com/2014/06/hypothesis-testing-of-temperature-trends.html

@author: brugger
"""

from climateUtils import getUAH, getGissLoti, getRSStlt
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import linregress
import statsmodels.api as sm

def autoCorr(x):
    rho = [np.sum(x*x), np.sum(x[1:]*x[:-1]), np.sum(x[2:]*x[:-2])]
    rho = rho / rho[0]
    return rho

def calcNu(x):
    rho = sm.tsa.acf(sm.tsa.detrend(x), nlags=2, fft=False)
    nu = 1 + 2*rho[1] / (1 - rho[2]/rho[1])
    return max(nu, 1)

d, t = getUAH()
#d, t = getGissLoti()
#d, t = getRSStlt()
mask = d >= 1979
d = d[mask]
t = t[mask]

result = linregress(d, t)
nu = calcNu(t)

print("{0:5.3f}, {1:5.3f}\n".format(10*result.slope, 10*2*np.sqrt(nu)*result.stderr))

for k in range(1,13):
    result = linregress(d[k::12], t[k::12])
    nu = calcNu(t[k::12])
    print("{0:d}, {1:5.3f}, {2:5.3f}".format(k, 10*result.slope, 10*2*np.sqrt(nu)*result.stderr))

# plt.plot(d, t, '-b')