mirror of
https://github.com/wassname/Volt.git
synced 2026-09-09 11:16:09 +08:00
tidy
This commit is contained in:
wassname
@@ -1,80 +0,0 @@
|
||||
import matplotlib.pyplot as plt
|
||||
import seaborn as sns
|
||||
import numpy as np
|
||||
import torch
|
||||
import pandas as pd
|
||||
import os
|
||||
import gpytorch
|
||||
import argparse
|
||||
import datetime
|
||||
|
||||
from botorch.models import SingleTaskGP
|
||||
from botorch.optim.fit import fit_gpytorch_torch
|
||||
from gpytorch.likelihoods import GaussianLikelihood
|
||||
from gpytorch.mlls import ExactMarginalLogLikelihood
|
||||
from gpytorch.means import ConstantMean, LinearMean
|
||||
from gpytorch.kernels import SpectralMixtureKernel, MaternKernel, RBFKernel, ScaleKernel
|
||||
from voltron.means import EWMAMean, DEWMAMean, TEWMAMean
|
||||
from voltron.train_utils import LearnGPCV, TrainVolModel, TrainVoltMagpieModel, TrainBasicModel
|
||||
from voltron.models import VoltMagpie
|
||||
from voltron.means import LogLinearMean
|
||||
|
||||
from voltron.rollout_utils import GeneratePrediction, Rollouts
|
||||
from voltron.data import make_ticker_list, DataGetter, GetStockHistory
|
||||
|
||||
|
||||
def BasicWindRollouts(train_x, train_y, test_x, kernel_name, mean_name='ewma', k=20,
|
||||
train_iters=600, nsample=1000):
|
||||
|
||||
|
||||
kernel_possibilities = {"sm": SpectralMixtureKernel,
|
||||
"matern": MaternKernel,
|
||||
"rbf": RBFKernel}
|
||||
kernel = kernel_possibilities[kernel_name.lower()]
|
||||
if kernel_name.lower() != "sm":
|
||||
kernel = ScaleKernel(kernel())
|
||||
else:
|
||||
kernel = kernel(num_mixtures=20)
|
||||
kernel.initialize_from_data_empspect(train_x, train_y.log())
|
||||
|
||||
model = SingleTaskGP(
|
||||
train_x.view(-1,1),
|
||||
train_y.log().reshape(-1, 1),
|
||||
covar_module=kernel,
|
||||
likelihood=GaussianLikelihood()
|
||||
)
|
||||
|
||||
mean_name = mean_name.lower()
|
||||
if mean_name == "loglinear":
|
||||
model.mean_module = LogLinearMean(1)
|
||||
model.mean_module.initialize_from_data(train_x, train_y.log())
|
||||
elif mean_name == 'linear':
|
||||
model.mean_module = LinearMean(1)
|
||||
elif mean_name == "constant":
|
||||
model.mean_module = ConstantMean()
|
||||
elif mean_name == "ewma":
|
||||
model.mean_module = EWMAMean(train_x, train_y.log(), k=k).to(train_x.device)
|
||||
elif mean_name == "dewma":
|
||||
model.mean_module = DEWMAMean(train_x, train_y.log(), k=k).to(train_x.device)
|
||||
elif mean_name == "tewma":
|
||||
model.mean_module = TEWMAMean(train_x, train_y.log(), k=k).to(train_x.device)
|
||||
|
||||
|
||||
model = model.to(train_x.device)
|
||||
mll = ExactMarginalLogLikelihood(model.likelihood, model)
|
||||
fit_gpytorch_torch(mll, options={'maxiter':train_iters, 'disp':False})
|
||||
|
||||
if mean_name in ["loglinear", "constant", 'linear']:
|
||||
save_samples = model.posterior(test_x).sample(torch.Size((nsample,
|
||||
))).squeeze(-1).cpu().detach()
|
||||
else:
|
||||
save_samples = Rollouts(
|
||||
train_x, train_y, test_x, model, nsample=nsample, method = "nonvol"
|
||||
).cpu().detach()
|
||||
|
||||
|
||||
torch.cuda.empty_cache()
|
||||
del model
|
||||
|
||||
|
||||
return save_samples
|
||||
@@ -1,178 +0,0 @@
|
||||
import numpy as np
|
||||
import torch
|
||||
import pandas as pd
|
||||
import gpytorch
|
||||
import argparse
|
||||
import datetime
|
||||
import warnings
|
||||
import copy
|
||||
import os
|
||||
from voltron.data import make_ticker_list, GetStockHistory
|
||||
import sys
|
||||
sys.path.append("../calibration")
|
||||
from LSTMUtils import SequenceDataset, LSTM, TrainLSTM, LSTMRollouts, NLL
|
||||
from torch.utils.data import DataLoader
|
||||
from voltron.train_utils import LearnGPCV, TrainVolModel, TrainVoltMagpieModel, TrainBasicModel
|
||||
from voltron.rollout_utils import Rollouts
|
||||
from BasicWind import BasicWindRollouts
|
||||
import pickle as pkl
|
||||
|
||||
def main(args):
|
||||
|
||||
stn_names, stn_lonlat, full_data = pkl.load(open("./wind_data.p", 'rb'))
|
||||
|
||||
use_cuda = False
|
||||
if torch.cuda.is_available():
|
||||
use_cuda = True
|
||||
|
||||
stn = args.stn_idx
|
||||
ntest = args.forecast_horizon
|
||||
ntrain = args.ntrain
|
||||
n_test_times = args.n_test_times
|
||||
ntime = full_data[0].shape[0]
|
||||
|
||||
test_idxs = torch.arange(ntrain, ntime-ntest,
|
||||
int((ntime-ntest-ntrain)/n_test_times))
|
||||
|
||||
stn_idxs = list(stn_names.keys())
|
||||
if args.kernel == 'volt':
|
||||
train_x = torch.arange(ntrain-1).float()/365
|
||||
else:
|
||||
train_x = torch.arange(ntrain).float()/365
|
||||
test_x = torch.arange(ntrain, ntrain + ntest).float()/365
|
||||
|
||||
if use_cuda:
|
||||
train_x, test_x = train_x.cuda(), test_x.cuda()
|
||||
|
||||
savepath = "./saved-outputs/stn" + str(stn) + "/"
|
||||
stn_data = full_data[stn]
|
||||
stn_data[stn_data == -99.0] = 0.
|
||||
if stn_data.mean() != 0:
|
||||
if not os.path.exists(savepath):
|
||||
os.mkdir(savepath)
|
||||
|
||||
for last_day in test_idxs:
|
||||
# try:
|
||||
raw_y = stn_data[last_day-ntrain:last_day] + 1
|
||||
train_y = torch.FloatTensor(raw_y)
|
||||
if use_cuda:
|
||||
train_y = train_y.cuda()
|
||||
|
||||
if args.kernel == 'volt':
|
||||
with gpytorch.settings.max_cholesky_size(2000):
|
||||
vol = LearnGPCV(train_x, train_y, train_iters=200,
|
||||
printing=False)
|
||||
vmod, vlh = TrainVolModel(train_x, vol,
|
||||
train_iters=500, printing=False)
|
||||
|
||||
if args.mean == 'constant':
|
||||
voltron, lh = TrainVoltMagpieModel(train_x, train_y[1:],
|
||||
vmod, vlh, vol,
|
||||
printing=False,
|
||||
train_iters=200, mean_func="constant")
|
||||
vmod.eval();
|
||||
voltron.eval();
|
||||
voltron.vol_model.eval();
|
||||
theta = 0.01
|
||||
# for theta in [0., 0.01, 0.025, 0.05, 0.1]:
|
||||
|
||||
temp_model = copy.deepcopy(voltron)
|
||||
with torch.no_grad():
|
||||
save_samples = Rollouts(train_x, train_y, test_x, temp_model,
|
||||
nsample=args.nsample, theta=theta)
|
||||
torch.save(save_samples, savepath + args.kernel + "_theta" + str(theta) +\
|
||||
"_" + str(last_day.item()) + ".pt")
|
||||
|
||||
del temp_model
|
||||
|
||||
else:
|
||||
for k in [400]:
|
||||
voltron, lh = TrainVoltMagpieModel(train_x, train_y[1:],
|
||||
vmod, vlh, vol,
|
||||
printing=False,
|
||||
train_iters=0, mean_func="ewma", k=k)
|
||||
vmod.eval();
|
||||
voltron.eval();
|
||||
voltron.vol_model.eval();
|
||||
for theta in [0.01]:
|
||||
temp_model = copy.deepcopy(voltron)
|
||||
with torch.no_grad():
|
||||
save_samples = Rollouts(train_x, train_y,
|
||||
test_x, temp_model,
|
||||
nsample=args.nsample, theta=theta)
|
||||
torch.save(save_samples, savepath + args.kernel + "_ema" + str(k) +\
|
||||
"_theta" + str(theta) +\
|
||||
"_" + str(last_day.item()) + ".pt")
|
||||
del temp_model
|
||||
del voltron, vmod, vol, vlh
|
||||
else:
|
||||
k=200
|
||||
rollouts = BasicWindRollouts(train_x, train_y, test_x,
|
||||
train_iters=args.train_epochs,
|
||||
kernel_name=args.kernel,
|
||||
mean_name=args.mean, k=k,
|
||||
nsample=200)
|
||||
|
||||
torch.save(rollouts, savepath + args.kernel + "_" +\
|
||||
args.mean + str(k) + "_" + str(last_day.item()) + ".pt")
|
||||
|
||||
print("stn ", stn, " idx ", last_day.item())
|
||||
|
||||
torch.cuda.empty_cache()
|
||||
# except:
|
||||
# print("### BROKEN stn", stn, " idx", last_day, " ###")
|
||||
if __name__ == "__main__":
|
||||
parser = argparse.ArgumentParser()
|
||||
parser.add_argument(
|
||||
"--stn_idx",
|
||||
type=int,
|
||||
default=0,
|
||||
)
|
||||
parser.add_argument(
|
||||
"--mean",
|
||||
type=str,
|
||||
default='constant',
|
||||
)
|
||||
parser.add_argument(
|
||||
"--n_test_times",
|
||||
type=int,
|
||||
default=10,
|
||||
)
|
||||
parser.add_argument(
|
||||
"--forecast_horizon",
|
||||
type=int,
|
||||
default=100,
|
||||
)
|
||||
parser.add_argument(
|
||||
"--kernel",
|
||||
type=str,
|
||||
default="matern",
|
||||
)
|
||||
parser.add_argument(
|
||||
"--ntrain",
|
||||
type=int,
|
||||
default=400,
|
||||
)
|
||||
parser.add_argument(
|
||||
"--nsample",
|
||||
type=int,
|
||||
default=1000,
|
||||
)
|
||||
parser.add_argument(
|
||||
"--printing",
|
||||
type=bool,
|
||||
default=False
|
||||
)
|
||||
parser.add_argument(
|
||||
"--train_epochs",
|
||||
type=int,
|
||||
default=500,
|
||||
)
|
||||
parser.add_argument(
|
||||
"--save",
|
||||
type=bool,
|
||||
default=False,
|
||||
)
|
||||
args = parser.parse_args()
|
||||
|
||||
main(args)
|
||||
@@ -1,15 +0,0 @@
|
||||
This directory contains the code needed to run Volt+Magpie on wind speed data taken from the [U.S. Climate Reference Network](https://www.ncei.noaa.gov/access/crn/).
|
||||
|
||||
|
||||
To source the data first walk through the `make_wind_dataset` notebook.
|
||||
|
||||
To generate forecasts for a station then run
|
||||
|
||||
```{bash}
|
||||
python GPGenerator.py
|
||||
--kernel={volt, sm, matern} ## kernel choice
|
||||
--stn_idx=0 ## station index in the dataset
|
||||
--mean={ewma, constant} ## mean choice
|
||||
--ntrain=400 ## training window
|
||||
--n_test_times=100 ## number of test time points
|
||||
```
|
||||
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
@@ -1,2 +0,0 @@
|
||||
robinhood_username="greg.w.benton@gmail.com"
|
||||
robinhood_password="Tho561mas!"
|
||||
@@ -1,12 +0,0 @@
|
||||
__version__ = 'alpha'
|
||||
from .kernels import BMKernel, VolatilityKernel
|
||||
from .models import BMGP, MultitaskBMGP
|
||||
from .train_utils import LearnGPCV
|
||||
from .option_utils import *
|
||||
try:
|
||||
from .robinhood_utils import GetStockData
|
||||
except:
|
||||
print("Warning no robinhood utils.")
|
||||
|
||||
from .rollout_utils import Rollouts, GeneratePrediction
|
||||
|
||||
@@ -1,41 +0,0 @@
|
||||
import numpy as np
|
||||
import torch
|
||||
import pandas as pd
|
||||
|
||||
|
||||
def GetTrainingData(SPY, date, N):
|
||||
idx = SPY[SPY["Date"] == date].index.item()
|
||||
return SPY['Close'].iloc[(idx-N):idx]
|
||||
|
||||
def GetTrueValue(SPY, date, strike):
|
||||
close_px = SPY['Close'][SPY["Date"] == date].item()
|
||||
return np.maximum(close_px-strike, 0)
|
||||
|
||||
def GetTradingDays(SPY, start, stop):
|
||||
start_idx = SPY[SPY["Date"] == start].index.item()
|
||||
stop_idx = SPY[SPY["Date"] == stop].index.item()
|
||||
return stop_idx-start_idx
|
||||
|
||||
def FindLastTradingDays(SPY, dates):
|
||||
last_days = []
|
||||
for date in dates:
|
||||
last_days.append(np.max(np.where(SPY.Date < date)[0]))
|
||||
|
||||
return np.array(SPY.Date[last_days])
|
||||
|
||||
def Pricer(mc_pxs, options, edays, true_pxs):
|
||||
logger = []
|
||||
for eday_idx, eday in enumerate(edays):
|
||||
eday = pd.Timestamp(eday)
|
||||
opts = options[options.expiration==pd.Timestamp(eday)]
|
||||
for idx, row in opts.iterrows():
|
||||
K = row.strike
|
||||
bid = row.bid
|
||||
ask = row.ask
|
||||
valuation = np.mean(np.maximum(mc_pxs[:, eday_idx].numpy() - K, 0))
|
||||
rtn = np.maximum(true_pxs[eday_idx] - K, 0)
|
||||
logger.append([eday, K, bid, ask, valuation, rtn.item()])
|
||||
|
||||
df = pd.DataFrame(logger)
|
||||
df.columns = ['Expiry', "Strike", "Bid", "Ask", "Voltron", "Return"]
|
||||
return df
|
||||
@@ -1,22 +0,0 @@
|
||||
import robin_stocks.robinhood as r
|
||||
import os
|
||||
import pandas as pd
|
||||
from dotenv import load_dotenv
|
||||
|
||||
def GetStockData(symbols, interval='day', span='5year'):
|
||||
"""
|
||||
just a wrapper for robin-stocks calls
|
||||
"""
|
||||
load_dotenv()
|
||||
username = os.getenv("robinhood_username")
|
||||
password = os.getenv("robinhood_password")
|
||||
r.login(username, password);
|
||||
|
||||
data = pd.DataFrame(r.stocks.get_stock_historicals(symbols, interval, span))
|
||||
data['date'] = pd.to_datetime(data['begins_at'], format='%Y-%m-%d').dt.date
|
||||
|
||||
ohlc = ['open_price', 'close_price', 'high_price', 'low_price']
|
||||
data[ohlc] = data[ohlc].astype("float")
|
||||
|
||||
return data[['date', 'symbol', 'open_price', 'close_price',
|
||||
'high_price', 'low_price']]
|
||||
@@ -1,115 +0,0 @@
|
||||
import torch
|
||||
import gpytorch
|
||||
from gpytorch.utils.cholesky import psd_safe_cholesky
|
||||
from gpytorch.utils.cholesky import psd_safe_cholesky
|
||||
|
||||
def GeneratePrediction(train_x, train_y, test_x, pred_vol, model, latent_mean=None, theta=0.5):
|
||||
vol = model.log_vol_path.exp()
|
||||
if train_x.ndim != test_x.ndim:
|
||||
test_x_for_stack = test_x.unsqueeze(0).repeat(train_x.shape[0], 1)
|
||||
else:
|
||||
test_x_for_stack = test_x
|
||||
if vol.ndim == 1:
|
||||
vol_for_stack = vol.unsqueeze(0).repeat(pred_vol.shape[0], 1)
|
||||
else:
|
||||
vol_for_stack = vol
|
||||
|
||||
full_x = torch.cat((train_x, test_x_for_stack),dim=-1)
|
||||
# print("vol stack = ", vol_for_stack.shape)
|
||||
# print("pred_vol = ", pred_vol.shape)
|
||||
full_vol = torch.cat((vol_for_stack, pred_vol),dim=-1)
|
||||
|
||||
test_x.repeat(2, test_x.numel())
|
||||
|
||||
idx_cut = train_x.shape[-1]
|
||||
|
||||
cov_mat = model.covar_module(full_x.unsqueeze(-1), full_vol.unsqueeze(-1)).evaluate()
|
||||
K_tr = cov_mat[..., :idx_cut, :idx_cut]
|
||||
K_tr_te = cov_mat[..., :idx_cut, idx_cut:]
|
||||
K_te = cov_mat[..., idx_cut:, idx_cut:]
|
||||
|
||||
train_mean = model.mean_module(train_x)
|
||||
train_diffs = train_y.unsqueeze(-1) - train_mean.unsqueeze(-1)
|
||||
# use psd cholesky if you must evaluate
|
||||
K_tr_chol = psd_safe_cholesky(K_tr, jitter=1e-4)
|
||||
pred_mean = K_tr_te.transpose(-1, -2).matmul(torch.cholesky_solve(train_diffs, K_tr_chol))
|
||||
# print(voltron.mean_module(test_x).detach().T.shape)
|
||||
# print(pred_mean.shape)
|
||||
pred_mean += model.mean_module(test_x).detach().T.unsqueeze(-1)
|
||||
|
||||
if latent_mean is not None:
|
||||
pred_mean -= theta * (pred_mean - latent_mean)
|
||||
|
||||
pred_cov = K_te - K_tr_te.transpose(-1, -2).matmul(torch.cholesky_solve(K_tr_te, K_tr_chol))
|
||||
|
||||
pred_cov_L = psd_safe_cholesky(pred_cov, jitter=1e-4)
|
||||
samples = torch.randn(*cov_mat.shape[:-2], test_x.shape[0], 1).to(test_x.device)
|
||||
samples = pred_cov_L @ samples
|
||||
|
||||
if pred_mean.ndim == 1:
|
||||
return samples + pred_mean.unsqueeze(-1)
|
||||
else:
|
||||
return (samples + pred_mean).squeeze(-1)
|
||||
|
||||
|
||||
|
||||
def Rollouts(train_x, train_y, test_x, model, nsample=50, method = "volt", theta=None,
|
||||
return_vol=False):
|
||||
if method != "volt":
|
||||
return nonvol_rollouts(train_x, train_y, test_x, model, nsample=nsample)
|
||||
if theta is None:
|
||||
latent_mean = None
|
||||
else:
|
||||
latent_mean = train_y.log().mean()
|
||||
ntest = test_x.numel()
|
||||
samples = torch.zeros(nsample, ntest)
|
||||
pred_vol = model.vol_model(test_x).sample(torch.Size((nsample, ))).exp()
|
||||
samples[:, 0] = GeneratePrediction(train_x, train_y,
|
||||
test_x[0].unsqueeze(0),
|
||||
pred_vol[:, 0].unsqueeze(1),
|
||||
model, latent_mean, theta).squeeze()
|
||||
train_stack_y = train_y.repeat(nsample, 1)
|
||||
train_stack_vol = model.log_vol_path.repeat(nsample, 1)
|
||||
|
||||
for idx in range(1, ntest):
|
||||
stack_y = torch.cat((train_stack_y,
|
||||
samples[:, :idx].to(train_stack_y.device)), -1)
|
||||
stack_vol = torch.cat((train_stack_vol,
|
||||
pred_vol[:, :idx].to(train_stack_vol.device).log()), -1)
|
||||
|
||||
rolling_x = torch.cat((train_x, test_x[:idx]))
|
||||
model.mean_module.train_y = stack_y
|
||||
model.mean_module.train_x = rolling_x
|
||||
|
||||
# train_x = rolling_x
|
||||
# train_y = stack_y
|
||||
model.log_vol_path = stack_vol
|
||||
samples[:, idx] = GeneratePrediction(rolling_x, stack_y,
|
||||
test_x[idx].unsqueeze(0),
|
||||
pred_vol[:, idx].unsqueeze(-1),
|
||||
model, latent_mean, theta).squeeze()
|
||||
if return_vol:
|
||||
return samples, pred_vol
|
||||
else:
|
||||
return samples
|
||||
|
||||
def nonvol_rollouts(train_x, train_y, test_x, model, nsample=50):
|
||||
ntest = test_x.numel()
|
||||
samples = torch.zeros(nsample, ntest)
|
||||
samples[:, 0] = model.posterior(test_x[0].unsqueeze(0)).sample(torch.Size((nsample,))).squeeze().squeeze()
|
||||
train_stack_y = train_y.repeat(nsample, 1)
|
||||
|
||||
for idx in range(1, ntest):
|
||||
stack_y = torch.cat((train_stack_y,
|
||||
samples[:, :idx].to(train_stack_y.device)), -1)
|
||||
rolling_x = torch.cat((train_x, test_x[:idx]))
|
||||
|
||||
model.mean_module.train_y = stack_y
|
||||
model.mean_module.train_x = rolling_x
|
||||
|
||||
model.train_inputs = (rolling_x.view(-1,1),)
|
||||
model.train_targets = stack_y
|
||||
model.train() # clear any caches that might have built up
|
||||
test_pt = test_x[idx].view(-1,1)
|
||||
samples[:, idx] = model.posterior(test_pt).sample().squeeze()
|
||||
return samples
|
||||
@@ -1,94 +0,0 @@
|
||||
import numpy as np
|
||||
import torch
|
||||
import gpytorch
|
||||
import sys
|
||||
|
||||
sys.path.append("../")
|
||||
from voltron.likelihoods import VolatilityGaussianLikelihood
|
||||
from voltron.models import SingleTaskVariationalGP
|
||||
from voltron.kernels import BMKernel, VolatilityKernel, FBMKernel
|
||||
from voltron.models import BMGP, BasicGP, Volt
|
||||
from voltron.means import LogLinearMean, EWMAMean, DEWMAMean, TEWMAMean, MeanRevertingEMAMean
|
||||
from gpytorch.kernels import ScaleKernel, RBFKernel, MaternKernel
|
||||
|
||||
|
||||
def LearnGPCV(train_x, train_y, train_iters=1000, printing=False, early_stopping=False, kernel = "bm"):
|
||||
dt = train_x[1]-train_x[0]
|
||||
scaled_returns = (train_y[1:] - train_y[:-1]) / (train_y[:-1]) / (dt**0.5)
|
||||
yy = scaled_returns
|
||||
|
||||
likelihood = VolatilityGaussianLikelihood(param="exp")
|
||||
# likelihood.raw_a.data -= 4.
|
||||
if kernel == "bm":
|
||||
covar_module = BMKernel()
|
||||
elif kernel == "fbm":
|
||||
covar_module = FBMKernel()
|
||||
model = SingleTaskVariationalGP(
|
||||
init_points=train_x.view(-1,1), likelihood=likelihood, use_piv_chol_init=False,
|
||||
mean_module = gpytorch.means.ConstantMean(), covar_module=covar_module,
|
||||
learn_inducing_locations=False, use_whitened_var_strat=False
|
||||
)
|
||||
model.initialize_variational_parameters(likelihood, train_x, y=yy)
|
||||
|
||||
model.train()
|
||||
likelihood.train()
|
||||
|
||||
# Use the adam optimizer
|
||||
optimizer = torch.optim.Adam([
|
||||
{"params": model.parameters()},
|
||||
# {"params": likelihood.parameters(), "lr": 0.1}
|
||||
], lr=0.01)
|
||||
|
||||
# "Loss" for GPs - the marginal log likelihood
|
||||
# num_data refers to the number of training datapoints
|
||||
mll = gpytorch.mlls.VariationalELBO(likelihood, model, yy.numel(), combine_terms = True)
|
||||
|
||||
print_every = 50
|
||||
for i in range(train_iters):
|
||||
# Zero backpropped gradients from previous iteration
|
||||
optimizer.zero_grad()
|
||||
# Get predictive output
|
||||
with gpytorch.settings.num_gauss_hermite_locs(75):
|
||||
output = model(train_x)
|
||||
# Calc loss and backprop gradients
|
||||
loss = -mll(output, yy)
|
||||
loss.backward()
|
||||
|
||||
if printing:
|
||||
if i % print_every == 0:
|
||||
print('Iter %d/%d - Loss: %.3f' % (i + 1, train_iters, loss.item()))
|
||||
optimizer.step()
|
||||
model.eval();
|
||||
likelihood.eval();
|
||||
predictive = model(train_x)
|
||||
pred_scale = likelihood(predictive, return_gaussian=False).scale.mean(0).detach()
|
||||
|
||||
return pred_scale
|
||||
|
||||
def TrainVolModel(train_x, vol_path, train_iters=1000, printing=False, kernel = "bm"):
|
||||
vol_lh = gpytorch.likelihoods.GaussianLikelihood().to(train_x.device)
|
||||
vol_lh.noise.data = torch.tensor([1e-2])
|
||||
vol_model = BMGP(train_x, vol_path.log(), vol_lh, kernel=kernel).to(train_x.device)
|
||||
# vol_model.covar_module.raw_vol.data = torch.tensor([-3.])
|
||||
|
||||
optimizer = torch.optim.Adam([
|
||||
{'params': vol_model.parameters()}, # Includes GaussianLikelihood parameters
|
||||
], lr=0.01)
|
||||
|
||||
# "Loss" for GPs - the marginal log likelihood
|
||||
mll = gpytorch.mlls.ExactMarginalLogLikelihood(vol_lh, vol_model)
|
||||
|
||||
print_every = 50
|
||||
for i in range(train_iters):
|
||||
# Zero gradients from previous iteration
|
||||
optimizer.zero_grad()
|
||||
# Output from model
|
||||
output = vol_model(train_x)
|
||||
# Calc loss and backprop gradients
|
||||
loss = -mll(output, vol_path.log())
|
||||
loss.backward()
|
||||
if printing:
|
||||
if i % print_every == 0:
|
||||
print('Iter %d/%d - Loss: %.3f' % (i + 1, train_iters, loss.item()))
|
||||
optimizer.step()
|
||||
return vol_model, vol_lh
|
||||
@@ -1,41 +0,0 @@
|
||||
import pandas as pd
|
||||
import yfinance as yf
|
||||
import datetime
|
||||
|
||||
|
||||
def make_ticker_list(file_name):
|
||||
tickers = open(file_name, 'r')
|
||||
tickers = [i.strip() for i in list(tickers)]
|
||||
return tickers
|
||||
|
||||
def make_price_files(tickers, start, end, fpath, printing):
|
||||
for i in tickers:
|
||||
history = yf.download(tickers=i,
|
||||
start=start,
|
||||
end=end,
|
||||
progress=False,
|
||||
)
|
||||
history.to_csv(fpath + str(i) + '.csv')
|
||||
if printing:
|
||||
print(str(i))
|
||||
|
||||
|
||||
def DataGetter(history = 500, fpath="../data/", printing=False, end_date=None,
|
||||
ticker_file="test_tickers.txt"):
|
||||
if end_date is None:
|
||||
end_date = datetime.date.today()
|
||||
else:
|
||||
end_date = datetime.datetime.strptime(end_date, "%Y-%m-%d").date()
|
||||
|
||||
start_date = end_date - datetime.timedelta(history)
|
||||
end_date = str(end_date)
|
||||
|
||||
tickers = make_ticker_list(fpath + ticker_file)
|
||||
make_price_files(tickers, start_date, end_date, fpath, printing)
|
||||
|
||||
def GetStockHistory(ticker, end_date=str(datetime.date.today()), history=500):
|
||||
end_date = datetime.datetime.strptime(end_date, "%Y-%m-%d").date()
|
||||
start_date = end_date - datetime.timedelta(history)
|
||||
return yf.download(tickers=ticker, start=start_date, end=end_date, progress=False)
|
||||
|
||||
|
||||
@@ -1 +0,0 @@
|
||||
from .MakeData import make_ticker_list, make_price_files, DataGetter, GetStockHistory
|
||||
@@ -1,10 +0,0 @@
|
||||
ADBE
|
||||
GOOG
|
||||
AMZN
|
||||
AMAT
|
||||
BRK-B
|
||||
DAL
|
||||
MCD
|
||||
NFLX
|
||||
PENN
|
||||
ZBRA
|
||||
@@ -1,16 +0,0 @@
|
||||
import torch
|
||||
from torch.nn.functional import softplus
|
||||
from gpytorch.kernels import Kernel
|
||||
|
||||
class BMKernel(Kernel):
|
||||
def __init__(self, vol=0., **kwargs):
|
||||
super(BMKernel, self).__init__(**kwargs)
|
||||
self.register_parameter(name='raw_vol',
|
||||
parameter=torch.nn.Parameter(vol*torch.ones(1)))
|
||||
|
||||
def forward(self, x1s, x2s, **kwargs):
|
||||
|
||||
X1, X2 = torch.meshgrid(x1s[:, 0], x2s[:, 0])
|
||||
# return self.raw_vol.exp() * torch.minimum(X1,X2)
|
||||
cov = self.raw_vol.exp() * torch.minimum(X1,X2)
|
||||
return cov
|
||||
@@ -1,37 +0,0 @@
|
||||
import torch
|
||||
from gpytorch.kernels import Kernel
|
||||
|
||||
def CumTrapz(y, x):
|
||||
dx = x[1] - x[0]
|
||||
wghts = dx * torch.ones_like(x)
|
||||
wghts[0] *= 0.5
|
||||
wghts[-1] *= 0.5
|
||||
return torch.cumsum(wghts * y, 0)
|
||||
|
||||
class VolatilityKernel(Kernel):
|
||||
has_lengthscale = False
|
||||
|
||||
def __init__(self, **kwargs):
|
||||
super().__init__(**kwargs)
|
||||
|
||||
def forward(self, x, vol_path, diag=False, **params):
|
||||
|
||||
last_dim_is_batch = params.get("last_dim_is_batch", False)
|
||||
if not last_dim_is_batch:
|
||||
vol_int = CumTrapz(vol_path.squeeze()**2, x.squeeze())
|
||||
else:
|
||||
x = x.unsqueeze(-1).repeat(x, vol_path.shape[-1])
|
||||
vol_int = CumTrapz(vol_path.pow(2.0), x)
|
||||
|
||||
idx = torch.arange(x.shape[0])
|
||||
idx1, idx2 = torch.meshgrid(idx, idx)
|
||||
idx = torch.minimum(idx1, idx2)
|
||||
res = vol_int[idx]
|
||||
|
||||
if vol_path.shape[-1] > 1:
|
||||
res = res.permute(2, 0, 1)
|
||||
|
||||
if diag:
|
||||
return torch.diagonal(res, dim1=-2, dim2=-1)
|
||||
else:
|
||||
return res
|
||||
@@ -1,2 +0,0 @@
|
||||
from .BMKernel import BMKernel
|
||||
from .VolKernel import VolatilityKernel
|
||||
@@ -1 +0,0 @@
|
||||
from .volatility_likelihood import VolatilityGaussianLikelihood
|
||||
@@ -1,61 +0,0 @@
|
||||
import torch
|
||||
|
||||
from torch.distributions import Normal
|
||||
from gpytorch.constraints import Positive, Interval
|
||||
from gpytorch.likelihoods import Likelihood, _OneDimensionalLikelihood
|
||||
|
||||
|
||||
class VolatilityGaussianLikelihood(_OneDimensionalLikelihood):
|
||||
def __init__(self, K=5, batch_shape=torch.Size(), param="cv", *args, **kwargs):
|
||||
"""
|
||||
parameterization of gaussian likelihood for volatility models like in
|
||||
wilson & ghahramani, copula processes, eq. 21.
|
||||
|
||||
we also consider the gp-exp parameterization
|
||||
"""
|
||||
|
||||
super().__init__()
|
||||
if param == "cv":
|
||||
self.raw_a = torch.nn.Parameter(torch.rand(*batch_shape, K, requires_grad=True))
|
||||
raw_b_init = 0.1 * torch.rand(*batch_shape, K)
|
||||
self.raw_b = torch.nn.Parameter(raw_b_init.detach_().requires_grad_())
|
||||
self.raw_c = torch.nn.Parameter(torch.rand(*batch_shape, K, requires_grad=True))
|
||||
|
||||
self.register_constraint("raw_a", Positive())
|
||||
self.register_constraint("raw_b", Interval(0.0, 3.0))
|
||||
self.register_constraint("raw_c", Interval(-3.0, 3.0))
|
||||
# elif param == "exp":
|
||||
# print("Using gp-exp parameterization.")
|
||||
self.param = param
|
||||
|
||||
@property
|
||||
def trans_a(self):
|
||||
return self.raw_a_constraint.transform(self.raw_a)
|
||||
|
||||
@property
|
||||
def trans_b(self):
|
||||
return self.raw_b_constraint.transform(self.raw_b)
|
||||
|
||||
@property
|
||||
def trans_c(self):
|
||||
return self.raw_c_constraint.transform(self.raw_c)
|
||||
|
||||
def forward(self, function_samples, *args, **kwargs):
|
||||
if self.param == "cv":
|
||||
transform = (
|
||||
(self.trans_b * function_samples.unsqueeze(-1) + self.trans_c).exp() + 1
|
||||
).log() * self.trans_a
|
||||
summed_transform = transform.sum(-1)
|
||||
else:
|
||||
summed_transform = function_samples.exp()
|
||||
return Normal(torch.zeros_like(summed_transform), summed_transform.clamp(min=1e-3))
|
||||
|
||||
def expected_log_prob(self, target, input, *params, **kwargs):
|
||||
res = super().expected_log_prob(target, input, *params, **kwargs)
|
||||
num_event_dim = len(input.event_shape)
|
||||
if num_event_dim > 1:
|
||||
res = res.sum(-1)
|
||||
return res
|
||||
|
||||
|
||||
# TODO: use a multitask Gaussian likelihood somehow in the multitask setting
|
||||
@@ -1,113 +0,0 @@
|
||||
import torch
|
||||
import gpytorch
|
||||
from gpytorch.means import Mean
|
||||
import numpy as np
|
||||
|
||||
def _EWMA(y, k):
|
||||
alpha = 2./(k + 1)
|
||||
conv = torch.nn.Conv1d(1, 1, kernel_size=k)
|
||||
wghts = alpha * (1-alpha)**(torch.arange(k-1, -1, -1))
|
||||
conv.weight.data = wghts.unsqueeze(0).unsqueeze(0)/wghts.sum()
|
||||
conv.bias.data = torch.zeros(1)
|
||||
|
||||
padded_px = torch.cat((y.squeeze()[0] * torch.ones(k),
|
||||
y.squeeze()))
|
||||
padded_px = padded_px.reshape(1, 1, -1)
|
||||
with torch.no_grad():
|
||||
ma = conv(padded_px).squeeze()
|
||||
return ma.type(torch.FloatTensor)
|
||||
|
||||
def EWMA(y, k):
|
||||
alpha = 2./(k + 1)
|
||||
conv = torch.nn.Conv1d(1, 1, kernel_size=k)
|
||||
wghts = alpha * (1-alpha)**(torch.arange(k-1, -1, -1))
|
||||
conv.weight.data = wghts.unsqueeze(0).unsqueeze(0)/wghts.sum()
|
||||
conv.bias.data = torch.zeros(1)
|
||||
|
||||
conv = conv.to(y.device)
|
||||
res = y[..., 0].unsqueeze(-1) * torch.ones(*y.shape[:-1], k).to(y.device)
|
||||
padded_px = torch.cat((res, y), dim=-1)
|
||||
batch_dim = y.shape[-2] if y.ndim > 1 else 1
|
||||
padded_px = padded_px.reshape(batch_dim, 1, -1)
|
||||
# print("padded_px shape = ", padded_px.shape)
|
||||
with torch.no_grad():
|
||||
ma = conv(padded_px).squeeze()
|
||||
|
||||
# print("ma shape = ", ma.shape)
|
||||
return ma.type(torch.FloatTensor)
|
||||
|
||||
class EWMAMean(Mean):
|
||||
def __init__(self, train_x, train_y, k=20):
|
||||
super().__init__()
|
||||
self.k = k
|
||||
self.train_x = train_x
|
||||
self.train_y = train_y
|
||||
|
||||
def forward(self, x):
|
||||
ewma = EWMA(self.train_y, self.k)
|
||||
if x.numel() == 1:
|
||||
res = ewma[..., -1].unsqueeze(0)
|
||||
return res.type(torch.FloatTensor).to(self.train_x.device)
|
||||
elif torch.equal(x.squeeze(), self.train_x.squeeze()):
|
||||
return ewma[..., :-1].type(torch.FloatTensor).to(self.train_x.device)
|
||||
else:
|
||||
return ewma.type(torch.FloatTensor).to(self.train_x.device)
|
||||
|
||||
|
||||
class HEWMAMean(Mean):
|
||||
def __init__(self, train_x, train_y, k=20):
|
||||
super().__init__()
|
||||
self.k = k
|
||||
self.train_x = train_x
|
||||
self.train_y = train_y
|
||||
|
||||
def forward(self, x):
|
||||
wma_k = EWMA(self.train_y, self.k)
|
||||
wma_k2 = EWMA(self.train_y, int(self.k/2))
|
||||
hma = EWMA(2*wma_k2[:-1] - wma_k[:-1], int(np.sqrt(self.k)))
|
||||
if torch.equal(x.squeeze(), self.train_x.squeeze()):
|
||||
return hma[:-1].type(torch.FloatTensor).to(self.train_x.device)
|
||||
else:
|
||||
return hma.type(torch.FloatTensor).to(self.train_x.device)
|
||||
|
||||
|
||||
class DEWMAMean(Mean):
|
||||
def __init__(self, train_x, train_y, k=20):
|
||||
super().__init__()
|
||||
self.k = k
|
||||
self.train_x = train_x
|
||||
self.train_y = train_y
|
||||
|
||||
def forward(self, x):
|
||||
ema = EWMA(self.train_y, self.k)#[..., :-1]
|
||||
ema_ema = EWMA(ema, self.k)[..., :-1]
|
||||
dema = 2*ema - ema_ema
|
||||
if x.numel() == 1:
|
||||
res = dema[..., -1].unsqueeze(0)
|
||||
return res.type(torch.FloatTensor).to(self.train_x.device)
|
||||
elif torch.equal(x.squeeze(), self.train_x.squeeze()):
|
||||
return dema[..., :-1].type(torch.FloatTensor).to(self.train_x.device)
|
||||
else:
|
||||
return dema.type(torch.FloatTensor).to(self.train_x.device)
|
||||
|
||||
|
||||
class TEWMAMean(Mean):
|
||||
def __init__(self, train_x, train_y, k=20):
|
||||
super().__init__()
|
||||
self.k = k
|
||||
self.alpha = 2./(self.k + 1)
|
||||
self.train_x = train_x
|
||||
self.train_y = train_y
|
||||
|
||||
def forward(self, x):
|
||||
ema = EWMA(self.train_y, self.k)
|
||||
ema_ema = EWMA(ema, self.k)[..., :-1]
|
||||
ema_ema_ema = EWMA(ema_ema, self.k)[..., :-1]
|
||||
tema = 3*ema - 3*ema_ema + ema_ema_ema
|
||||
if x.numel() == 1:
|
||||
res = tema[..., -1].unsqueeze(0)
|
||||
return res.type(torch.FloatTensor).to(self.train_x.device)
|
||||
elif torch.equal(x.squeeze(), self.train_x.squeeze()):
|
||||
return tema[..., :-1].type(torch.FloatTensor).to(self.train_x.device)
|
||||
else:
|
||||
return tema.type(torch.FloatTensor).to(self.train_x.device)
|
||||
@@ -1,3 +0,0 @@
|
||||
from .loglinear_mean import LogLinearMean
|
||||
from .mulidentity_mean import MulIdentityMean
|
||||
from .EWMA import EWMAMean, DEWMAMean, TEWMAMean
|
||||
@@ -1,21 +0,0 @@
|
||||
import torch
|
||||
|
||||
from gpytorch.means import LinearMean
|
||||
|
||||
class LogLinearMean(LinearMean):
|
||||
def __init__(self, input_size, batch_shape=None, bias=True):
|
||||
if batch_shape is None:
|
||||
batch_shape = torch.Size()
|
||||
|
||||
super().__init__(input_size=input_size, batch_shape=batch_shape, bias=bias)
|
||||
|
||||
def initialize_from_data(self, x, y):
|
||||
with torch.no_grad():
|
||||
# assume y is on log scale
|
||||
self.bias.data = y.exp().mean(-1)
|
||||
# is there anything we should do for the mean term?
|
||||
|
||||
def forward(self, x):
|
||||
linear_term = super().forward(x)
|
||||
# to prevent linear stuff
|
||||
return linear_term.clamp(min=1e-6).log()
|
||||
@@ -1,53 +0,0 @@
|
||||
import math
|
||||
import torch
|
||||
import gpytorch
|
||||
import numpy as np
|
||||
from voltron.means import EWMAMean, DEWMAMean, TEWMAMean
|
||||
from botorch.models import SingleTaskGP
|
||||
from botorch.optim.fit import fit_gpytorch_torch
|
||||
from voltron.rollout_utils import nonvol_rollouts
|
||||
|
||||
class BasicGP():
|
||||
def __init__(self, train_x, train_y, kernel="matern", mean='constant',
|
||||
k=400, num_mixtures=10):
|
||||
# super(BasicGP, self).__init__(train_x, train_y, likelihood)
|
||||
if mean.lower() == 'constant':
|
||||
mean_module = gpytorch.means.ConstantMean().to(train_x.device)
|
||||
elif mean.lower() == 'ewma':
|
||||
mean_module = EWMAMean(train_x, train_y, k).to(train_x.device)
|
||||
elif mean.lower() == 'dewma':
|
||||
mean_module = DEWMAMean(train_x, train_y, k).to(train_x.device)
|
||||
elif mean.lower() == 'tewma':
|
||||
mean_module = TEWMAMean(train_x, train_y, k).to(train_x.device)
|
||||
else:
|
||||
print("ERROR: Mean not implemented")
|
||||
|
||||
if kernel.lower() == 'matern':
|
||||
covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.MaternKernel())
|
||||
elif kernel.lower() in ['sm', 'spectralmixture', 'spectral']:
|
||||
covar_module = gpytorch.kernels.SpectralMixtureKernel(num_mixtures=num_mixtures)
|
||||
covar_module.initialize_from_data(train_x, train_y)
|
||||
elif kernel.lower() == 'rbf':
|
||||
covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel())
|
||||
else:
|
||||
print("ERROR: Kernel not implemented")
|
||||
|
||||
self.model = SingleTaskGP(train_x.view(-1, 1), train_y.reshape(-1, 1),
|
||||
covar_module=covar_module,
|
||||
likelihood=gpytorch.likelihoods.GaussianLikelihood())
|
||||
self.model.mean_module = mean_module
|
||||
|
||||
def Train(self, train_iters=400, display=False):
|
||||
mll = gpytorch.mlls.ExactMarginalLogLikelihood(self.model.likelihood, self.model)
|
||||
fit_gpytorch_torch(mll, options={'maxiter':train_iters, 'disp':display})
|
||||
|
||||
def Forecast(self, test_x, nsample=100):
|
||||
if not isinstance(self.model.mean_module, (EWMAMean, DEWMAMean, TEWMAMean)):
|
||||
|
||||
samples = self.model.posterior(test_x).sample(torch.Size((nsample, )))
|
||||
|
||||
else:
|
||||
samples = nonvol_rollouts(self.model.train_inputs[0].squeeze(),
|
||||
self.model.train_targets.squeeze(),
|
||||
test_x, self.model, nsample)
|
||||
return samples.squeeze()
|
||||
@@ -1,113 +0,0 @@
|
||||
import numpy as np
|
||||
import torch
|
||||
from torch import nn
|
||||
from torch.utils.data import DataLoader, Dataset
|
||||
from torch.autograd import Variable
|
||||
|
||||
|
||||
class SequenceDataset(Dataset):
|
||||
def __init__(self, data, sequence_length=5):
|
||||
self.sequence_length = sequence_length
|
||||
self.X = data.float()
|
||||
|
||||
def __len__(self):
|
||||
return self.X.shape[0]-1
|
||||
|
||||
def __getitem__(self, i):
|
||||
if i >= self.sequence_length - 1:
|
||||
i_start = i - self.sequence_length + 1
|
||||
x = self.X[i_start:(i + 1)]
|
||||
else:
|
||||
padding = self.X[0].repeat(self.sequence_length - i - 1, 1).squeeze(-1)
|
||||
x = self.X[0:(i + 1)]
|
||||
x = torch.cat((padding, x), 0)
|
||||
|
||||
return x.unsqueeze(0), self.X[i+1]
|
||||
|
||||
class LSTM(nn.Module):
|
||||
def __init__(self, train_x, train_y, seq_len, hidden_size,
|
||||
num_layers, batch_size=128):
|
||||
super(LSTM, self).__init__()
|
||||
|
||||
self.train_x = train_x
|
||||
self.train_y = train_y
|
||||
|
||||
self.norm_y = (train_y - train_y.mean())/train_y.std()
|
||||
|
||||
self.dset = SequenceDataset(self.norm_y, sequence_length=seq_len)
|
||||
self.trainloader = DataLoader(self.dset, batch_size=batch_size,
|
||||
shuffle=True)
|
||||
|
||||
self.num_classes = 1
|
||||
self.num_layers = num_layers
|
||||
self.input_size = seq_len
|
||||
self.hidden_size = hidden_size
|
||||
|
||||
self.lstm = nn.LSTM(input_size=seq_len, hidden_size=hidden_size,
|
||||
num_layers=num_layers, batch_first=True) #lstm
|
||||
self.fc_1 = nn.Linear(hidden_size, 128) #fully connected 1
|
||||
self.fc = nn.Linear(128, 2) #fully connected last layer
|
||||
|
||||
self.relu = nn.ReLU()
|
||||
self.softplus = nn.Softplus()
|
||||
|
||||
def forward(self,x):
|
||||
h_0 = Variable(torch.zeros(self.num_layers, x.size(0), self.hidden_size)).to(x.device) #hidden state
|
||||
c_0 = Variable(torch.zeros(self.num_layers, x.size(0), self.hidden_size)).to(x.device) #internal state
|
||||
# Propagate input through LSTM
|
||||
output, (hn, cn) = self.lstm(x, (h_0, c_0)) #lstm with input, hidden, and internal state
|
||||
|
||||
hn = hn[self.num_layers-1]
|
||||
hn = hn.view(-1, self.hidden_size) #reshaping the data for Dense layer next
|
||||
out = self.relu(hn)
|
||||
out = self.fc_1(out) #first Dense
|
||||
out = self.relu(out) #relu
|
||||
out = self.fc(out) #Final Output
|
||||
|
||||
output = torch.zeros_like(out)
|
||||
output[:, 0] = out[:, 0]
|
||||
output[:, 1] = self.softplus(out[:, 1])
|
||||
return output
|
||||
|
||||
def Loss(self, targets, outputs):
|
||||
dist = torch.distributions.Normal(outputs[:, 0], outputs[:, 1])
|
||||
return -dist.log_prob(targets).sum()
|
||||
|
||||
def Train(self, epochs, display=False):
|
||||
optimizer = torch.optim.Adam(self.parameters(), lr=0.01)
|
||||
num_batches = len(self.trainloader)
|
||||
total_loss = 0
|
||||
self.train()
|
||||
for epoch in range(epochs):
|
||||
for X, y in self.trainloader:
|
||||
X = X.to(self.train_x.device)
|
||||
y = y.to(self.train_x.device)
|
||||
output = self(X)
|
||||
loss = self.Loss(y, output)
|
||||
|
||||
optimizer.zero_grad()
|
||||
loss.backward()
|
||||
optimizer.step()
|
||||
|
||||
total_loss += loss.item()
|
||||
if display:
|
||||
if epoch%50 == 0:
|
||||
avg_loss = total_loss / num_batches
|
||||
print(f"Train loss: {avg_loss}, Epoch: {epoch}")
|
||||
|
||||
|
||||
def Forecast(self, test_x, nsample=50):
|
||||
rollout_len = test_x.shape[0]
|
||||
xin, xout = self.dset[len(self.dset)-1]
|
||||
xx = torch.cat((xin[0, 1:], xout.unsqueeze(0)))
|
||||
xx = xx.repeat(nsample, 1).unsqueeze(1)
|
||||
xx = xx.to(self.train_x.device)
|
||||
roll_pxs = torch.zeros(nsample, rollout_len)
|
||||
with torch.no_grad():
|
||||
for idx in range(rollout_len):
|
||||
out = self(xx)
|
||||
smpl = torch.normal(out[:, 0], out[:, 1])
|
||||
roll_pxs[:, idx] = smpl
|
||||
xx = torch.cat((xx[..., 1:], smpl.unsqueeze(-1).unsqueeze(-1)), -1)
|
||||
return roll_pxs * self.train_y.std() + self.train_y.mean()
|
||||
|
||||
@@ -1,162 +0,0 @@
|
||||
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
|
||||
@@ -1,6 +0,0 @@
|
||||
from .BMGP import BMGP, MultitaskBMGP
|
||||
from .multi_task_variational_gp import MultitaskVariationalGP
|
||||
from .single_task_variational_gp import SingleTaskVariationalGP
|
||||
from .BasicGPModels import BasicGP
|
||||
from .Volt import Volt
|
||||
from .LSTM import LSTM
|
||||
@@ -1,265 +0,0 @@
|
||||
from typing import Union
|
||||
from copy import deepcopy
|
||||
|
||||
import torch
|
||||
import functools
|
||||
|
||||
from botorch.models.gpytorch import GPyTorchModel
|
||||
from botorch.models import SingleTaskGP
|
||||
from botorch.posteriors import GPyTorchPosterior
|
||||
|
||||
from gpytorch import lazify
|
||||
from gpytorch.distributions import MultivariateNormal
|
||||
from gpytorch.lazy import (
|
||||
CholLazyTensor,
|
||||
TriangularLazyTensor,
|
||||
)
|
||||
from gpytorch.likelihoods import GaussianLikelihood
|
||||
from gpytorch.likelihoods import FixedNoiseGaussianLikelihood as FNGaussianLikelihood
|
||||
from gpytorch.likelihoods.gaussian_likelihood import _GaussianLikelihoodBase
|
||||
from gpytorch.means import ConstantMean
|
||||
from gpytorch.models import ApproximateGP
|
||||
from gpytorch.kernels import ScaleKernel, RBFKernel, InducingPointKernel
|
||||
from gpytorch.utils.errors import NotPSDError
|
||||
from gpytorch.utils.memoize import cached, add_to_cache, clear_cache_hook
|
||||
from gpytorch.variational import (
|
||||
CholeskyVariationalDistribution,
|
||||
UnwhitenedVariationalStrategy,
|
||||
VariationalStrategy,
|
||||
)
|
||||
|
||||
# from ..utils import pivoted_cholesky_init
|
||||
|
||||
# copied from wjmaddox/volatilitygp
|
||||
|
||||
# def _update_caches(m, *args, **kwargs):
|
||||
# if hasattr(m, "_memoize_cache"):
|
||||
# for key, item in m._memoize_cache.items():
|
||||
# if type(item) is not tuple and type(item) is not MultivariateNormal:
|
||||
# if len(args) is 0:
|
||||
# new_lc = item.to(torch.empty(0, **kwargs))
|
||||
# else:
|
||||
# new_lc = item.to(*args)
|
||||
# m._memoize_cache[key] = new_lc
|
||||
# if type(item) is TriangularLazyTensor:
|
||||
# m._memoize_cache[key] = m._memoize_cache[key].double()
|
||||
# elif type(item) is MultivariateNormal:
|
||||
# if len(args) is 0:
|
||||
# new_lc = item.lazy_covariance_matrix.to(torch.empty(0, **kwargs))
|
||||
# else:
|
||||
# new_lc = item.lazy_covariance_matrix.to(*args)
|
||||
# m._memoize_cache[key] = MultivariateNormal(
|
||||
# item.mean.to(*args, **kwargs), new_lc
|
||||
# )
|
||||
# else:
|
||||
# m._memoize_cache[key] = (x.to(*args, **kwargs) for x in item)
|
||||
|
||||
|
||||
# def _add_cache_hook(tsr, pred_strat):
|
||||
# if tsr.grad_fn is not None:
|
||||
# wrapper = functools.partial(clear_cache_hook, pred_strat)
|
||||
# functools.update_wrapper(wrapper, clear_cache_hook)
|
||||
# tsr.grad_fn.register_hook(wrapper)
|
||||
# return tsr
|
||||
|
||||
|
||||
class _SingleTaskVariationalGP(ApproximateGP):
|
||||
def __init__(
|
||||
self,
|
||||
init_points: torch.Tensor = None,
|
||||
likelihood=None,
|
||||
learn_inducing_locations=True,
|
||||
covar_module=None,
|
||||
mean_module=None,
|
||||
use_piv_chol_init=True,
|
||||
num_inducing=None,
|
||||
use_whitened_var_strat=True,
|
||||
init_targets=None,
|
||||
train_inputs=None,
|
||||
train_targets=None,
|
||||
):
|
||||
|
||||
if covar_module is None:
|
||||
covar_module = ScaleKernel(RBFKernel())
|
||||
|
||||
# if use_piv_chol_init:
|
||||
# if num_inducing is None:
|
||||
# num_inducing = int(init_points.shape[-2] / 2)
|
||||
|
||||
# if num_inducing < init_points.shape[-2]:
|
||||
# covar_module = covar_module.to(init_points)
|
||||
|
||||
# covariance = covar_module(init_points)
|
||||
# if init_targets is not None and init_targets.shape[-1] == 1:
|
||||
# init_targets = init_targets.squeeze(-1)
|
||||
# if likelihood is not None and not isinstance(
|
||||
# likelihood, GaussianLikelihood
|
||||
# ):
|
||||
# _ = likelihood.newton_iteration(
|
||||
# init_points, init_targets, model=None, covar=covariance
|
||||
# )
|
||||
# if likelihood.has_diag_hessian:
|
||||
# hessian_sqrt = likelihood.expected_hessian().sqrt()
|
||||
# else:
|
||||
# hessian_sqrt = (
|
||||
# lazify(likelihood.expected_hessian())
|
||||
# .root_decomposition()
|
||||
# .root
|
||||
# )
|
||||
# covariance = hessian_sqrt.matmul(covariance).matmul(
|
||||
# hessian_sqrt.transpose(-1, -2)
|
||||
# )
|
||||
# inducing_points = pivoted_cholesky_init(
|
||||
# init_points, covariance.evaluate(), num_inducing
|
||||
# )
|
||||
# else:
|
||||
# inducing_points = init_points.detach().clone()
|
||||
# else:
|
||||
inducing_points = init_points.detach().clone()
|
||||
|
||||
variational_distribution = CholeskyVariationalDistribution(
|
||||
inducing_points.shape[-2]
|
||||
)
|
||||
if use_whitened_var_strat:
|
||||
variational_strategy = VariationalStrategy(
|
||||
self,
|
||||
inducing_points,
|
||||
variational_distribution,
|
||||
learn_inducing_locations=learn_inducing_locations,
|
||||
)
|
||||
else:
|
||||
variational_strategy = UnwhitenedVariationalStrategy(
|
||||
self,
|
||||
inducing_points,
|
||||
variational_distribution,
|
||||
learn_inducing_locations=learn_inducing_locations,
|
||||
)
|
||||
super(_SingleTaskVariationalGP, self).__init__(variational_strategy)
|
||||
self.mean_module = ConstantMean() if mean_module is None else mean_module
|
||||
self.mean_module.to(init_points)
|
||||
self.covar_module = covar_module
|
||||
|
||||
self.likelihood = GaussianLikelihood() if likelihood is None else likelihood
|
||||
self.likelihood.to(init_points)
|
||||
self.train_inputs = [train_inputs] if train_inputs is not None else [init_points]
|
||||
self.train_targets = train_targets if train_targets is not None else init_targets
|
||||
|
||||
self.condition_into_exact = True
|
||||
|
||||
self.to(init_points)
|
||||
|
||||
def forward(self, x):
|
||||
mean_x = self.mean_module(x)
|
||||
covar_x = self.covar_module(x)
|
||||
latent_pred = MultivariateNormal(mean_x, covar_x)
|
||||
return latent_pred
|
||||
|
||||
# may actually want to keep this one in the future
|
||||
# def to(self, *args, **kwargs):
|
||||
# _update_caches(self, *args, **kwargs)
|
||||
# self.variational_strategy = self.variational_strategy.to(*args, **kwargs)
|
||||
# _update_caches(self.variational_strategy, *args, **kwargs)
|
||||
# return super().to(*args, **kwargs)
|
||||
|
||||
|
||||
class SingleTaskVariationalGP(_SingleTaskVariationalGP, GPyTorchModel):
|
||||
def __init__(
|
||||
self,
|
||||
init_points=None,
|
||||
likelihood=None,
|
||||
learn_inducing_locations=True,
|
||||
covar_module=None,
|
||||
mean_module=None,
|
||||
use_piv_chol_init=True,
|
||||
num_inducing=None,
|
||||
use_whitened_var_strat=True,
|
||||
init_targets=None,
|
||||
train_inputs=None,
|
||||
train_targets=None,
|
||||
outcome_transform=None,
|
||||
input_transform=None,
|
||||
):
|
||||
if outcome_transform is not None:
|
||||
is_gaussian_likelihood = (
|
||||
isinstance(likelihood, GaussianLikelihood) or likelihood is None
|
||||
)
|
||||
if train_targets is not None and is_gaussian_likelihood:
|
||||
if train_targets.ndim == 1:
|
||||
train_targets = train_targets.unsqueeze(-1)
|
||||
train_targets, _ = outcome_transform(train_targets)
|
||||
|
||||
if init_targets is not None and is_gaussian_likelihood:
|
||||
init_targets, _ = outcome_transform(init_targets)
|
||||
init_targets = init_targets.squeeze(-1)
|
||||
|
||||
if train_targets is not None:
|
||||
train_targets = train_targets.squeeze(-1)
|
||||
|
||||
# unlike in the exact gp case we need to use the input transform to pre-define the inducing pts
|
||||
if input_transform is not None:
|
||||
if init_points is not None:
|
||||
init_points = input_transform(init_points)
|
||||
|
||||
_SingleTaskVariationalGP.__init__(
|
||||
self,
|
||||
init_points=init_points,
|
||||
likelihood=likelihood,
|
||||
learn_inducing_locations=learn_inducing_locations,
|
||||
covar_module=covar_module,
|
||||
mean_module=mean_module,
|
||||
use_piv_chol_init=use_piv_chol_init,
|
||||
num_inducing=num_inducing,
|
||||
use_whitened_var_strat=use_whitened_var_strat,
|
||||
init_targets=init_targets,
|
||||
train_inputs=train_inputs,
|
||||
train_targets=train_targets,
|
||||
)
|
||||
|
||||
if input_transform is not None:
|
||||
self.input_transform = input_transform.to(
|
||||
self.variational_strategy.inducing_points
|
||||
)
|
||||
|
||||
if outcome_transform is not None:
|
||||
self.outcome_transform = outcome_transform.to(
|
||||
self.variational_strategy.inducing_points
|
||||
)
|
||||
|
||||
def forward(self, x):
|
||||
x = self.transform_inputs(x)
|
||||
return super().forward(x)
|
||||
|
||||
@property
|
||||
def num_outputs(self) -> int:
|
||||
# we should only be able to have one output without a multitask variational strategy here
|
||||
return 1
|
||||
|
||||
# might be useful in the future though
|
||||
# def posterior(
|
||||
# self,
|
||||
# X: torch.Tensor,
|
||||
# observation_noise: Union[bool, torch.Tensor] = False,
|
||||
# **kwargs,
|
||||
# ):
|
||||
# if observation_noise and not isinstance(self.likelihood, _GaussianLikelihoodBase):
|
||||
# noiseless_posterior = super().posterior(
|
||||
# X=X, observation_noise=False, **kwargs
|
||||
# )
|
||||
# noiseless_mvn = noiseless_posterior.mvn
|
||||
# neg_hessian_f = self.likelihood.neg_hessian_f(noiseless_mvn.mean)
|
||||
# try:
|
||||
# likelihood_cov = neg_hessian_f.inverse()
|
||||
# except:
|
||||
# eye_like_hessian = torch.eye(
|
||||
# neg_hessian_f.shape[-2],
|
||||
# device=neg_hessian_f.device,
|
||||
# dtype=neg_hessian_f.dtype,
|
||||
# )
|
||||
# likelihood_cov = lazify(neg_hessian_f).inv_matmul(eye_like_hessian)
|
||||
|
||||
# noisy_mvn = type(noiseless_mvn)(
|
||||
# noiseless_mvn.mean, noiseless_mvn.lazy_covariance_matrix + likelihood_cov
|
||||
# )
|
||||
# return GPyTorchPosterior(mvn=noisy_mvn)
|
||||
|
||||
# return super().posterior(X=X, observation_noise=observation_noise, **kwargs)
|
||||
Reference in New Issue
Block a user