<jupyter_start><jupyter_text><jupyter_code># SEt environment Variable
import os

#os.environ['QUANDI_API_KEY']=''
#!export | grep QUANDI
!pip install quandl
import warnings
import itertools
import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
import os

import quandl

from statsmodels.tsa.stattools import adfuller
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf


from sklearn.metrics import mean_squared_error
from statsmodels.tsa.arima_model import ARIMA
from statsmodels.tsa.arima_model import ARIMAResults
from math import sqrt
from scipy.stats import boxcox


plt.style.use('fivethirtyeight')

%matplotlib inline


# Data from quandi API. returns a numpy array
quandl.ApiConfig.api_key = os.environ.get('QUANDI_API_KEY')

data = quandl.get('NASDAQOMX/XQC', start_date='2017-04-30', end_date='2019-04-29', returns="numpy")
data[:20]
datadf = pd.DataFrame(data)  

print(datadf.describe())
datadf.head()
datadf['Index Value'].plot()
plt.show()

# create a series

series = pd.Series(datadf['Index Value'].values, index=datadf['Trade Date'])

series.plot()
plt.show()

#group plot with years

groups = series.groupby(pd.TimeGrouper('A'))

years = pd.DataFrame()

plt.figure(figsize = (10,10))

n_groups = len(groups)
i = 1

for name, group in groups:
  plt.subplot((n_groups * 100) + 10 + i)
  i +=1
  plt.plot(group)
  
plt.show()
print(years)
print(groups)
# histograme

plt.figure(figsize = (10,10))
plt.subplot(211)
series.hist()
plt.subplot(212)
series.plot(kind='kde')
plt.show()
# helper function


# create a differenced series
def difference(dataset, interval=1):
	diff = list()
	for i in range(interval, len(dataset)):
		value = dataset[i] - dataset[i - interval]
		diff.append(value)
	return diff



# invert differenced value
def inverse_difference(history, yhat, interval=1):
	return yhat + history[-interval]


 
# evaluate an ARIMA model for a given order (p,d,q) and return RMSE
def evaluate_arima_model(X, arima_order):
	# prepare training dataset
	X = X.astype('float32')
	train_size = int(len(X) * 0.50)
	train, test = X[0:train_size], X[train_size:]
	history = [x for x in train]
	# make predictions
	predictions = list()
	for t in range(len(test)):
		# difference data
		weeks_in_year = 52
		diff = difference(history, weeks_in_year)
		model = ARIMA(diff, order=arima_order)
		model_fit = model.fit(trend='nc', disp=0)
		yhat = model_fit.forecast()[0]
		yhat = inverse_difference(history, yhat, weeks_in_year)
		predictions.append(yhat)
		history.append(test[t])
	# calculate out of sample error
	mse = mean_squared_error(test, predictions)
	rmse = sqrt(mse)
	return rmse





# evaluate combinations of p, d and q values for an ARIMA model
def evaluate_models(dataset, p_values, d_values, q_values):
	dataset = dataset.astype('float32')
	best_score, best_cfg = float("inf"), None
	for p in p_values:
		for d in d_values:
			for q in q_values:
				order = (p,d,q)
				try:
					mse = evaluate_arima_model(dataset, order)
					if mse < best_score:
						best_score, best_cfg = mse, order
					print('ARIMA%s RMSE=%.3f' % (order,mse))
				except:
					continue
	print('Best ARIMA%s RMSE=%.3f' % (best_cfg, best_score))

# create train and text split

# CREATE TRAIN-TEST DATA


split_point = len(series) - 100

dataset, validation = series[0:split_point], series[split_point:]
print('Dataset %d, Validation %d' % (len(dataset), len(validation)))
dataset.to_csv('dataset.csv')
validation.to_csv('validation.csv')
# prepare data
X = series.values
X = X.astype('float32')
train_size = int(len(X) * 0.50)

train, test = X[0:train_size], X[train_size:]
# walk-forward validation
history = [x for x in train]
predictions = list()
for i in range(len(test)):
	# difference data
	weeks_in_year = 52
	diff = difference(history, weeks_in_year)
	# predict
	model = ARIMA(diff, order=(2,1,0))
	model_fit = model.fit(trend='nc', disp=0)
	yhat = model_fit.forecast()[0]
	yhat = inverse_difference(history, yhat, weeks_in_year)
	predictions.append(yhat)
	# observation
	obs = test[i]
	history.append(obs)
	print('>Predicted=%.3f, Expected=%3.f' % (yhat, obs))
# report performance
mse = mean_squared_error(test, predictions)
rmse = sqrt(mse)
print('RMSE: %.3f' % rmse)
## CODE FROM BROWNLEE FOR SAVING MODEL

# monkey patch around bug in ARIMA class
def __getnewargs__(self):
	return ((self.endog),(self.k_lags, self.k_diff, self.k_ma))

ARIMA.__getnewargs__ = __getnewargs__


# prepare data
X = series.values
X = X.astype('float32')
# difference data
weeks_in_year = 52
diff = difference(X, weeks_in_year)
# fit model
model = ARIMA(diff, order=(2,0,0))
model_fit = model.fit(trend='nc', disp=0)
# bias constant, could be calculated from in-sample mean residual - CHANGED TO ZERO
bias = 0
# save model
model_fit.save('model.pkl')
np.save('model_bias.npy', [bias])

# CODE FOR MAKING PREDICTION USING SAVED MODEL - LOOKING MANY STEPS AHEAD

 
# prepare datasets

X = dataset.values.astype('float32')
history = [x for x in X]
weeks_in_year = 52
validation = pd.Series.from_csv('validation.csv', header=0)
y = validation.values.astype('float32')

# load model
model_fit = ARIMAResults.load('model.pkl')
bias = np.load('model_bias.npy')

# make first prediction
predictions = list()
yhat = float(model_fit.forecast()[0])
yhat = bias + inverse_difference(history, yhat, weeks_in_year)
predictions.append(yhat)
history.append(y[0])
print('>Predicted=%.3f, Expected=%3.f' % (yhat, y[0]))

plt.figure(figsize=(10,10))
# rolling forecasts
for i in range(1, len(y)):
	# difference data
	weeks_in_year = 52
	diff = difference(history, weeks_in_year)
	# predict
	model = ARIMA(diff, order=(2,0,0))
	model_fit = model.fit(trend='nc', disp=0)
	yhat = model_fit.forecast()[0]
	yhat = bias + inverse_difference(history, yhat, weeks_in_year)
	predictions.append(yhat)
	# observation
	obs = y[i]
	history.append(obs)
	print('>Predicted=%.3f, Expected=%3.f' % (yhat, obs))
# report performance
mse = mean_squared_error(y, predictions)
rmse = sqrt(mse)
print('RMSE: %.3f' % rmse)
plt.plot(y)
plt.plot(predictions, color='red')
plt.show()

plt.figure(figsize=(10,10))
plot_acf(series, ax=plt.gca())
plt.show()

#pyplot.subplot(212)

plt.figure(figsize=(10,10))
plot_pacf(series, ax=plt.gca())
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1360329009807023, 'fraction_numerical': 0.019297690604239165, 'mean_word_length': 5.2480237154150196, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '593531', 'license': 'MIT', 'n_tokens_mistral': 2544, 'n_tokens_neox': 2353, 'n_words': 570, 'path': 'time_series_api_quandi_stock_analysis.ipynb', 'repo_name': 'toraaglobal/CaseStudies', 'size': '6322'}