Volt/voltron/models/Volt.py

import torch
from torch.nn.functional import softplus
import gpytorch
from gpytorch.kernels import Kernel
from gpytorch.means import ConstantMean
from gpytorch.utils.cholesky import psd_safe_cholesky

from voltron.models.BMGP import BMGP, MultitaskBMGP
from voltron.kernels import VolatilityKernel

# import sys
# sys.path.append("../means/")
from voltron.means import EWMAMean, DEWMAMean, TEWMAMean
from voltron.train_utils import LearnGPCV, TrainVolModel
from voltron.rollout_utils import Rollouts

class Volt(gpytorch.models.ExactGP):
    def __init__(self, train_x, log_data, mean='constant',
                 vol_path=None, k=25):
        
        # WE ASSUME IN THE BATCHED CASE THAT
        # TRAIN_X: N 
        # TRAIN_Y: T X N
        # VOL_PATH: T X N

        likelihood = gpytorch.likelihoods.GaussianLikelihood()
        
        super(Volt, self).__init__(train_x[1:], log_data[1:], likelihood)
        
        if log_data.ndim > 1:
            batch_shape = log_data.shape[:-1]
        else:
            batch_shape = torch.Size()
            
        if mean.lower() == 'constant':
            mean_module = gpytorch.means.ConstantMean().to(train_x.device)
        elif mean.lower() == 'ewma':
            mean_module = EWMAMean(train_x[1:], log_data[1:], k).to(train_x.device)
        elif mean.lower() == 'dewma':
            mean_module = DEWMAMean(train_x[1:], log_data[1:], k).to(train_x.device)
        elif mean.lower() == 'tewma':
            mean_module = TEWMAMean(train_x[1:], log_data[1:], k).to(train_x.device)
        else:
            print("ERROR: Mean not implemented")
            
        self.mean_module = mean_module.to(train_x.device)
        self.covar_module = VolatilityKernel().to(train_x.device)
        
        # but we store a T X N X 1 copy of train_x to maintain consistency w/ 
        # gpytorch
        if log_data.ndim > 1:
            self.train_x = train_x.unsqueeze(0).repeat(*batch_shape, 1)
        else:
            self.train_x = train_x
        self.train_y = log_data
        
        if vol_path is None:
            self.log_vol_path = -1 * torch.ones(train_x.shape[0]-1)
        else:
            self.log_vol_path = vol_path.log()
        
        self.train_cov = self.covar_module(self.train_x.unsqueeze(-1), self.log_vol_path.exp().unsqueeze(-1)).detach()
        
        if batch_shape == torch.Size():
            self.vol_lh = gpytorch.likelihoods.GaussianLikelihood()
            self.vol_model = BMGP(train_x, self.log_vol_path, self.vol_lh) 
        else:
            self.vol_lh = gpytorch.likelihoods.MultitaskGaussianLikelihood(num_tasks=batch_shape[0])
            self.vol_lh.noise = 1e-3
            # we want the vol path GP to be N x T shaped and train_x to be N shaped
            self.vol_model = MultitaskBMGP(train_x, self.log_vol_path.t(), self.vol_lh)
        
    def UpdateVolPath(self, vol_path):
        self.log_vol_path = vol_path.log()
        self.train_cov = self.covar_module(self.train_inputs[0], self.log_vol_path.exp())
        return
        
    def VolMLL(self):
        vol_mll = gpytorch.mlls.ExactMarginalLogLikelihood(self.vol_lh, self.vol_model)
        outputs = self.vol_model(self.train_x)
        return vol_mll(outputs, self.log_vol_path)
        
    def forward(self, x):
        mean_x = self.mean_module(x)
        if torch.equal(x, self.train_inputs[0]):
            covar_x = self.train_cov
#             print("TRAIN COV")
        else:
            covar_x = self.covar_module(x, self.log_vol_path.exp())
#             print("NOT TRAIN COV")
            
#         print(covar_x.evaluate().shape)
        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)
    
    def Train(self, gpcv_iters=400, vol_mod_iters=1000, data_mod_iters=400, display=False):
        x = self.train_x.squeeze()
        data = self.train_y.exp()

        ##############################
        ## Train GPCV and Vol Model ##
        ##############################
        vol = LearnGPCV(x[1:], data, gpcv_iters, printing=display)
        vmod, vlh = TrainVolModel(x[1:], vol, vol_mod_iters, printing=display)
        
        self.UpdateVolPath(vol)
        ######################
        ## Train Data Model ##
        ######################
        if isinstance(self.mean_module, (EWMAMean, DEWMAMean, TEWMAMean)):
            grad_flags = [True, False, False, False]
        else:
            grad_flags = [True, True, False, False, False]
            
        
        self.likelihood.raw_noise.data = torch.tensor([1e-5]).to(x.device)
        self.vol_lh = vlh.to(x.device)
        self.vol_model = vmod.to(x.device)
    
        for idx, p in enumerate(self.parameters()):
            p.requires_grad = grad_flags[idx]
            
        self.train();
        self.vol_lh.train();
        self.vol_model.train();
        
        
        optimizer = torch.optim.Adam([
            {'params': self.parameters()},  # Includes GaussianLikelihood parameters
        ], lr=0.1)
        mll = gpytorch.mlls.ExactMarginalLogLikelihood(self.likelihood, self)
        
        print_every = 50
        for i in range(data_mod_iters):
            # Zero gradients from previous iteration
            optimizer.zero_grad()
            # Output from model
            output = self(x[1:])
#             print(output)
#             print(data.log().shape)
            # Calc loss and backprop gradients
            loss = -mll(output, data.log()[1:])
            loss.backward()
            if display:
                if i % print_every == 0:
                    print('Iter %d/%d - Loss: %.3f' % (i + 1, data_mod_iters, loss.item()))
            optimizer.step()
            
        
    def Forecast(self, test_x, nsample=50, return_vol=False, mean_revert=False, theta=0.05):
        self.vol_model.eval();
        self.eval();
        latent_mean = None
        if mean_revert:
            latent_mean = self.train_targets.squeeze().mean()
        samples = Rollouts(self.train_inputs[0].squeeze(),
                           self.train_targets.squeeze(), 
                           test_x, self, 
                            nsample=nsample,
                          return_vol=return_vol,
                          latent_mean=latent_mean, theta=theta)
        
        return samples