#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat Jul 24 08:14:11 2021

Attempt to implement trend analysis described at
https://hpklima.blogspot.com/2014/06/confidence-intervals-around-temperature.html?m=0

@author: mark
"""
import numpy as np
from climateUtils import getUAH

def ols(x, y):
    b = ( (np.sum(x * y) - np.mean(y) * np.sum(x)) /
         (np.sum(x**2) - np.mean(x) * np.sum(x)) )
    a = np.mean(y) - b*np.mean(x)
    return b, a

def variance(x):
    return np.mean( (x - np.mean(x))**2 )

def covariance(x,y):
    return np.mean( (x - np.mean(x)) * (y - np.mean(y)) )
    
def ols_c(x, y):
    b = covariance(x,y) / variance(x)
    a = np.mean(y) - b*np.mean(x)
    return b, a

def sse(x, y, b, a):
    return np.sum( (y - (a + b*x))**2 )

def se_regression(x, y, b, a):
    return sse(x, y, b, a) / (len(x) - 2)

def sxx(x):
    return np.sum( (x - np.mean(x))**2 )

def sigma_sq_slope(x, y, b, a):
    return se_regression(x, y, b, a) / sxx(x)

def sigma_sq_slope_C(x, y, b, a):
    return nu(y) * sigma_sq_slope(x, y, b, a)

# Autocorrelationn at offsets 0 to n
def rho(y, n=2):
    rho = np.zeros(n+1)
    scale = np.correlate(y,y)[0]
    rho[0] = 1
    for lag in range(1,n+1):
        rho[lag] = np.correlate(y[lag:],y[:-lag])[0] / scale
    return rho
    
def nu(y):
    r = rho(y, n=2)
    return 1 + (2*r[1]) / (1 - r[2]/r[1])
        
if __name__ == "__main__":
    d, t = getUAH()
    b, a = ols_c(d, t)
    Se = np.sqrt( se_regression(d, t, b, a))
    sigma_sq = sigma_sq_slope(d, t, b, a)
    print((b, a), sigma_sq)