#!/usr/bin/env python
# encoding: utf-8
"""
Title: Cash flow models

Stochastic cash flow models for predicting the stumpage prices for Finnish
timber assortments using a monthly time series of 1986-2008. Adopted from
matlab code by Kari Hyytiäinen 2008.

Copyright (c) 2009 Simosol Oy. Distributed under the GNU General Public License 2.0.
"""

import os
import math
import copy
import numpy

from base_prices import PRICES, PRICECLASSES
from base_prices_2009 import PRICES_2009, PRICECLASSES_2009
from base_biomass_prices import BM_PRICES, BM_PRICECLASSES


# Timber parameters, estimated Saeed Bayazidi from monthly timber prices (same
# data as in base_prices). The three sets of parameters are associated with
# three alternative price models (or procedures), named 2, 5 and 8 (functions
# _eq_2, _eq_5 and _eq_8)
ETHA2 = numpy.array([
    0.1460,
    0.2045,
    0.0402,
    0.2615,
    0.1544,
    0.1035,
    0.3146], dtype=float)
XBAR2 = numpy.array([
    56.37,
    48.59,
    53.72,
    17.28,
    24.97,
    17.61,
    37.02], dtype=float)
NEWSIGMA2 = numpy.array([
    0.04658,
    0.05145,
    0.06727,
    0.05312,
    0.08753,
    0.04942,
    0.05352], dtype=float)

ALPHA5 = numpy.array([0,0,0,0,0,0,0], dtype=float)
BETA5 = numpy.array([
    0.0021393,
    0.0023501,
    0.0022107,
    0.0020146,
    0.0037554,
    0.0024931,
    0.0055998], dtype=float)
SIGMA5 = numpy.array([
    0.0042102**0.5,
    0.0046643**0.5,
    0.0041568**0.5,
    0.0040810**0.5,
    0.0077516**0.5,
    0.0050742**0.5,
    0.011239**0.5], dtype=float)

ALPHA8 = numpy.array([
    0.0084368,
    0.0066304,
    0.0220265,
    0.000824791,
    -0.0045725,
    -0.000642836,
    0.0089528], dtype=float)
BETA8 = numpy.array([0,0,0,0,0,0,0], dtype=float)
SIGMA8 = numpy.array([
    0.0042114**0.5,
    0.0046657**0.5,
    0.0041596**0.5,
    0.004082**0.5,
    0.0077575**0.5,
    0.0050754**0.5,
    0.0112468**0.5], dtype=float)


# Timber parameters, estimated Saeed Bayazidi from monthly timber prices (same
# data as in base_prices_2009). The parameters are associated with procedure
# OU (Ornstein-Uhlenbeck) and used in function _eq_Ornstein_Uhlenbeck
ETHA_OU = numpy.array([
    0.016019632,
    0.011082181,
    0.006080398,
    0.007519813,
    0.00682831,
    0.0189168], dtype=float)
PBAR_OU = numpy.array([
    49.19449158,
    46.65414378,
    42.66793045,
    14.05466388,
    21.56931488,
    15.77907114], dtype=float)
SIGMA_OU = numpy.array([
    0.00554103 ** 0.5,
    0.00816552 ** 0.5,
    0.00491644 ** 0.5,
    0.0057944 ** 0.5,
    0.00567155 ** 0.5,
    0.014743 ** 0.5], dtype=float)


# biomass parameters from Saeed
BM_ETHA2 = numpy.array([
    0.1544,
    0.1035,
    0.1035], dtype=float)
BM_XBAR2 = numpy.array([
    24.97,
    17.61,
    17.61], dtype=float)
BM_NEWSIGMA2 = numpy.array([
    0.08753,
    0.04942,
    0.04942], dtype=float)

BM_ALPHA5 = numpy.array([0,0,0], dtype=float)
BM_BETA5 = numpy.array([
    0.0037554,
    0.0024931,
    0.0024931], dtype=float)
BM_SIGMA5 = numpy.array([
    0.0077516**0.5,
    0.0050742**0.5,
    0.0050742**0.5], dtype=float)

BM_ALPHA8 = numpy.array([
    -0.0045725,
    -0.000642836,
    -0.000642836], dtype=float)
BM_BETA8 = numpy.array([0,0,0], dtype=float)
BM_SIGMA8 = numpy.array([
    0.0077575**0.5,
    0.0050754**0.5,
    0.0050754**0.5], dtype=float)


def _load_prices(prices):
    """Load the pickled price info, calculate incremental change in the
    stumpage prices and normalize the changes
    """
    ind = numpy.arange(0, prices.shape[0]-12, 1)
    # calculate the incremental proportional change of stumpage prices in one
    # year intervals
    delta = numpy.log10(prices[ind+12,:]) - numpy.log10(prices[ind,:])
    mean = numpy.mean(delta, axis=0)
    std = numpy.std(delta, axis=0)
    # return the actual prices and the standardized price changes
    return prices, (delta - mean) / std

def _eq_2(P, X, n, BM):
    """Stochastic stumpage price equation 2
    """
    if BM:
        for i in xrange(0, n-1, 1):
            P[i+1] = P[i] + BM_ETHA2 * (BM_XBAR2 - P[i]) + BM_NEWSIGMA2 * P[i] * X[i]
    else:
        for i in xrange(0, n-1, 1):
            P[i+1] = P[i] + ETHA2 * (XBAR2 - P[i]) + NEWSIGMA2 * P[i] * X[i]
    return P

def _eq_5(P, X, n, BM):
    """Stochastic stumpage price equation 5
    """
    if BM:
        for i in xrange(0, n-1, 1):
            P[i+1] = P[i] + (BM_ALPHA5 + BM_BETA5 * P[i]) + BM_SIGMA5 * P[i] * X[i]
    else:
        for i in xrange(0, n-1, 1):
            P[i+1] = P[i] + (ALPHA5 + BETA5 * P[i]) + SIGMA5 * P[i] * X[i]
    return P

def _eq_8(P, X, n, BM):
    """Stochastic stumpage price equation 8
    """
    if BM:
        for i in xrange(0, n-1, 1):
            P[i+1] = P[i] + (BM_ALPHA8 + BM_BETA8 * P[i]) + BM_SIGMA8 * P[i] * X[i]
    else:
        for i in xrange(0, n-1, 1):
            P[i+1] = P[i] + (ALPHA8 + BETA8 * P[i]) + SIGMA8 * P[i] * X[i]
    return P

def _eq_Ornstein_Uhlenbeck(P, X, n):
    """Stochastic stumpage price equation, Ornstein-Uhlenbeck method
    
    dpt = etha(pbar - pt)dt + sigma * dzt
    """
    for i in xrange(0, n-1, 1):
        P[i+1] = P[i] + ETHA_OU * (PBAR_OU - P[i]) + SIGMA_OU * X[i]
    return P

def _simulate(prices, delta, n, proc, BM, priceclasses):
    """Generate correlated multivariate random variables from N(0,1)
    distribution and generate a set of random prices for n years with the
    given procedure
    """
    mu = numpy.zeros(priceclasses, dtype=float)
    sigma = numpy.cov(delta, rowvar=0)
    X = numpy.random.multivariate_normal(mu, sigma, n-1)
    P = numpy.zeros((n, priceclasses), dtype=float)
    # the initial values of each assortment should be the means from the
    # assortments in the price data
    P[0,:] = prices.mean(axis=0)
    if proc == 2:
        P = _eq_2(P, X, n, BM)
    elif proc == 5:
        P = _eq_5(P, X, n, BM)
    elif proc == 8:
        P = _eq_8(P, X, n, BM)
    elif proc == 'OU':
        P = _eq_Ornstein_Uhlenbeck(P, X, n)
    return P


# CONSTANTS
# cash flow table classifiers for timber (species, assortment)
PINELOG = (1,1)
SPRUCELOG = (2,1)
BIRCHLOG1 = (3,1)
BIRCHLOG2 = (4,1)
PINEPULP = (1,2)
SPRUCEPULP = (2,2)
BIRCHPULP1 = (3,2)
BIRCHPULP2 = (4,2)
# cash flow table classifiers for biomass (species, bm_assortment)
PINESTUMP_ROOTS = (1,1)
PINEBRANCH = (1,2)
PINEFOLIAGE = (1,3)
SPRUCESTUMP_ROOTS = (2,1)
SPRUCEBRANCH = (2,2)
SPRUCEFOLIAGE = (2,3)
BIRCHSTUMP_ROOTS = (3,1)
BIRCHBRANCH = (3,2)
BIRCHFOLIAGE = (3,3)
BIRCH2STUMP_ROOTS = (4,1)
BIRCH2BRANCH = (4,2)
BIRCH2FOLIAGE = (4,3)


def random_prices(table, date, iteration, variables, parameters, errors):
    """Simulate random prices using GMR method
    """
    ret = True
    
    # the initial year is now fixed to 2008
    n = date.year - 2008 + 1
    BM = False
    BM = 'bm_assortment' in table.classifiers
    
    if BM:
        #prices of bm_assortments
        index = (0,1)
        
        # load the prices and compute standardized delta (increment), simulate
        # prices for the next n years
        # price array contents: stump_roots, branch, foliage
        prices, delta_st = _load_prices(BM_PRICES)
        P = _simulate(prices, delta_st, n, 2, BM, BM_PRICECLASSES)
    
        # store the prices of the last year
        table.set_values(index, PINESTUMP_ROOTS, P[-1,0])
        table.set_values(index, SPRUCESTUMP_ROOTS, P[-1,0])
        table.set_values(index, BIRCHSTUMP_ROOTS, P[-1,0])
        table.set_values(index, BIRCH2STUMP_ROOTS, P[-1,0])
        table.set_values(index, PINEBRANCH, P[-1,1])
        table.set_values(index, SPRUCEBRANCH, P[-1,1])
        table.set_values(index, BIRCHBRANCH, P[-1,1])
        table.set_values(index, BIRCH2BRANCH, P[-1,1])
        table.set_values(index, PINEFOLIAGE, P[-1,2])
        table.set_values(index, SPRUCEFOLIAGE, P[-1,2])
        table.set_values(index, BIRCHFOLIAGE, P[-1,2])
        table.set_values(index, BIRCH2FOLIAGE, P[-1,2])
    
    else:
        #prices of timber assortments
        index = (1,2)
        
        # load the prices and compute standardized delta (increment), simulate
        # prices for the next n years
        # price array contents: mat, kut, kot, mak, kuk, kok, all
        prices, delta_st = _load_prices(PRICES)
        P = _simulate(prices, delta_st, n, 2, BM, PRICECLASSES)
    
        # store the prices of the last year
        table.set_values(index, PINELOG, P[-1,0])
        table.set_values(index, SPRUCELOG, P[-1,1])
        table.set_values(index, BIRCHLOG1, P[-1,2])
        table.set_values(index, BIRCHLOG2, P[-1,2])
        table.set_values(index, PINEPULP, P[-1,3])
        table.set_values(index, SPRUCEPULP, P[-1,4])
        table.set_values(index, BIRCHPULP1, P[-1,5])
        table.set_values(index, BIRCHPULP2, P[-1,5])

    return ret

def random_prices_Ornstein_Uhlenbeck(table, date, iteration, variables,
                                     parameters, errors):
    """Simulate random prices using the Ornstein-Uhlenbeck method
    """
    ret = True
    
    # the initial year is now fixed to 2008
    n = date.year - 2009 + 1
    BM = False
    BM = 'bm_assortment' in table.classifiers
    
    if BM:
        #prices of bm_assortments
        index = (0,1)
        
        # load the prices and compute standardized delta (increment), simulate
        # prices for the next n years
        prices, delta_st = _load_prices(BM_PRICES)
        P = _simulate(prices, delta_st, n, 2, BM, BM_PRICECLASSES)
    
        # store the prices of the last year
        table.set_values(index, PINESTUMP_ROOTS, P[-1,0])
        table.set_values(index, SPRUCESTUMP_ROOTS, P[-1,0])
        table.set_values(index, BIRCHSTUMP_ROOTS, P[-1,0])
        table.set_values(index, BIRCH2STUMP_ROOTS, P[-1,0])
        table.set_values(index, PINEBRANCH, P[-1,1])
        table.set_values(index, SPRUCEBRANCH, P[-1,1])
        table.set_values(index, BIRCHBRANCH, P[-1,1])
        table.set_values(index, BIRCH2BRANCH, P[-1,1])
        table.set_values(index, PINEFOLIAGE, P[-1,2])
        table.set_values(index, SPRUCEFOLIAGE, P[-1,2])
        table.set_values(index, BIRCHFOLIAGE, P[-1,2])
        table.set_values(index, BIRCH2FOLIAGE, P[-1,2])
    
    else:
        #prices of timber assortments
        index = (1,2)
        
        # load the prices and compute standardized delta (increment), simulate
        # prices for the next n years
        prices, delta_st = _load_prices(PRICES_2009)
        P = _simulate(prices, delta_st, n, 'OU', BM, PRICECLASSES_2009)
    
        # store the prices of the last year
        table.set_values(index, PINELOG, P[-1,0])
        table.set_values(index, SPRUCELOG, P[-1,1])
        table.set_values(index, BIRCHLOG1, P[-1,2])
        table.set_values(index, BIRCHLOG2, P[-1,2])
        table.set_values(index, PINEPULP, P[-1,3])
        table.set_values(index, SPRUCEPULP, P[-1,4])
        table.set_values(index, BIRCHPULP1, P[-1,5])
        table.set_values(index, BIRCHPULP2, P[-1,5])

    return ret


def main():
    pass

if __name__ == '__main__':
    main()
