MATLAB VWAP, Part II: Backtest

After reading some Coding Horror (this specifically), I decided I might as well begin sharing more results. Albeit not terribly useful from a trading perspective--I'll explain shortly--, here are some results utilizing a small database of GS minute values (approximately just short of two weeks).  The backtest uses an extension of the VWAP code from yesterday following the simple rules of long when the 100 minute VWAP > 500 minute VWAP and short when 100 min VWAP < 500 min VWAP. A pretty simple, standard, elementary test.

Here are the results (with image):

Return: +13.36% (rounded)

PLEASE DO NOT TRADE ON THESE RESULTS. 

As I mentioned above, this should not be used from a trading perspective for several reasons:

1) It is run on far too few data points. Perhaps it just happened to work perfectly for this interval. I have no clue if this holds. (If you would like to donate historical data, please let me know!)

2) This is done during earnings, and, as you would see if you ran this for other stocks, could be hit or miss, especially with overnight holding (and, again, especially through earnings).

3) There are plenty of stocks this does not work for. Some with similar sized returns, although negative.

4) My data could be flawed, my routine should be suspect and you should always perform tests read on the internet for yourself!

5) I am sure there are 10 more reasons.

Now, here is the "backtest" code in function form. It's for MATLAB and requires that you have some data in  a MySQL database (very easy if you use my IB Historical Data Downloader/Importer):

function sum_log_change = vwap_backtest(symbol_string)

clear close_mat
clear volume_mat
clear vwap
clear vwap_2
clear buy_sell_matrix
clear trade
clear trade_prices
clear sum_log_change

vwap_length = 100;
vwap_length_2 = 500;

format long

symbol_string

%connect to the database CHANGE PASSWORD FROM ***** TO WHATEVER YOUR PASSWORD IS!!!
conn = database('stocks','root','****','com.mysql.jdbc.Driver','jdbc:mysql://localhost:3306/stocks');

%query the symbols in the database
str_query_close = ['SELECT close, volume FROM stock_prices_minute WHERE symbol = ''',symbol_string,  ''''];
close_cell = fetch(conn, str_query_close);

close_mat = cell2mat(close_cell(:,1));
volume_mat = cell2mat(close_cell(:,2));

close_times_volume = close_mat .* volume_mat;


for i = vwap_length:length(close_times_volume)
    volume_total = 0;
    volume_price_sum = 0;
    for j = 0:(vwap_length - 1)
        volume_total = volume_total + volume_mat(i-j);
        volume_price_sum = volume_price_sum + volume_mat(i-j)*close_mat(i-j);
    end
    vwap(i-(vwap_length - 1)) = volume_price_sum / volume_total;
end


for i = vwap_length_2:length(close_times_volume)
    volume_total = 0;
    volume_price_sum = 0;
    for j = 0:(vwap_length_2 - 1)
        volume_total = volume_total + volume_mat(i-j);
        volume_price_sum = volume_price_sum + volume_mat(i-j)*close_mat(i-j);
    end
    vwap_2(i-(vwap_length_2 - 1)) = volume_price_sum / volume_total;
end


%do comparing
for i = 1:(length(close_times_volume) - (vwap_length_2 - 1))
    if vwap_2(i) < vwap(i+(vwap_length_2 - vwap_length))
        buy_sell_matrix(i) = 1;
    else
        buy_sell_matrix(i) = -1;
    end
end

buy_sell_matrix;

%calculate trade prices
for i = 1:(length(close_times_volume) - (vwap_length_2 - 1))
    if i == 1
        trades(i) = close_mat(i + vwap_length_2);
    else
        if buy_sell_matrix(i) ~= buy_sell_matrix(i-1)
            trades(i) = close_mat(i + vwap_length_2);
        elseif i == (length(close_times_volume)-(vwap_length_2 - 1))
            trades(i) = close_mat(length(close_mat));
        else
            trades(i) = 0;
        end
    end
end

trades;

j = 1;
%matrix of prices trades took place at
for i = 1:length(trades)
    if trades(i) ~= 0
        trade_prices(j) = trades(i);
        j = j + 1;
    end
end


trade_prices;

log_change = 0;
sum_log_change = 0;
%now you can calculate lognormal return
for i = 1:(length(trade_prices) - 1)
    if buy_sell_matrix(1) == 1
        if mod(i,2) == 0
            log_change = -(log(trade_prices(i+1)) - log(trade_prices(i)));
        else
            log_change = log(trade_prices(i+1)) - log(trade_prices(i));
        end
        sum_log_change = sum_log_change + log_change;
    else
        if mod(i,2)==0
            log_change = log(trade_prices(i+1)) - log(trade_prices(i));
        else
            log_change = -(log(trade_prices(i+1)) - log(trade_prices(i)));
        end
        sum_log_change = sum_log_change + log_change;
    end
end

sum_log_change;
            
        
    

x = [1:length(close_mat)];
x2 = [(vwap_length):length(close_mat)];
x3 = [(vwap_length_2):length(close_mat)];
plot(x, close_mat)
hold on
plot(x2, vwap, '-r')
hold on
plot(x3, vwap_2, '-g')

Link to previous article: MATLAB VWAP (Part I)

MATLAB VWAP

Here's a simple MATLAB script to extract data from a MySQL Database (here, specifically one that was created using my IB historical data extractor) and plot close data along with two VWAPs (volume weighted average prices) of your desired length. It's pretty simple but could be expanded upon to backtest VWAP strategies using the created VWAP matrices. (Also, I think I fixed the way source code is displayed on this site...should now be easier to copy and paste).

Update: Here's part II.

clear close_mat
clear volume_mat
clear vwap
clear vwap_2

vwap_length = 100;
vwap_length_2 = 500;

format long

%connect to the database
conn = database('stocks','root','***','com.mysql.jdbc.Driver','jdbc:mysql://localhost:3306/stocks');

%query the symbols in the database
str_query_close = 'SELECT close, volume FROM stock_prices_minute WHERE symbol = ''GOOG''';
close_cell = fetch(conn, str_query_close);

close_mat = cell2mat(close_cell(:,1))
volume_mat = cell2mat(close_cell(:,2));

close_times_volume = close_mat .* volume_mat;


for i = vwap_length:length(close_times_volume)
    volume_total = 0;
    volume_price_sum = 0;
    for j = 0:(vwap_length - 1)
        volume_total = volume_total + volume_mat(i-j);
        volume_price_sum = volume_price_sum + volume_mat(i-j)*close_mat(i-j);
    end
    vwap(i-(vwap_length - 1)) = volume_price_sum / volume_total;
end

for i = vwap_length_2:length(close_times_volume)
    volume_total = 0;
    volume_price_sum = 0;
    for j = 0:(vwap_length_2 - 1)
        volume_total = volume_total + volume_mat(i-j);
        volume_price_sum = volume_price_sum + volume_mat(i-j)*close_mat(i-j);
    end
    vwap_2(i-(vwap_length_2 - 1)) = volume_price_sum / volume_total;
end


x = [1:length(close_mat)];
x2 = [(vwap_length):length(close_mat)];
x3 = [(vwap_length_2):length(close_mat)];
plot(x, close_mat)
hold on
plot(x2, vwap, '-r')
hold on
plot(x3, vwap_2, '-g')

Downloading Options Data From Yahoo!

I wanted to do some volatility surfaces using Yahoo! options data. Unfortunately they don't make it as easy to grab options data as it is to grab stock data.

Here's a very nice solution. (To prevent loss, source code follows; although all credit goes to link!)

function DataOut = Get_Yahoo_Options_Data(symbolid)
%Get_Yahoo_Options_Data get Option Chain Data from Yahoo
% Get Options Chain Data from Yahoo
% DataOut = Get_Yahoo_Options_Data(symbol)
% Inputs: Symbol name as a character String
% Output:  A structure with the following fields
%       data : A 1xN cell where N is the number of Expiries available
%       ExpDates : A 1xN cell array of Expiry Dates
%       Calls  : A 1xN cell array of Call Option data for each expiry
%       Puts  : A 1xN cell array of Put Option data
%       CPHeaders : Headers for the calls and puts option data
%       Headers: Headers for the data
%       FullOptionData : A combined cell array of DataOut.data
%       Last : Last Price
% Example:
%           DataOut = Get_Yahoo_Options_Data('LVS');
% (c)tradingwithmatlab.blogspot.com
DataOut = struct;
% Construct and read the URL from Yahoo Finance Website
urlText = urlread(['http://finance.yahoo.com/q/os?s=' symbolid]);
% Try getting the Table Data from URL Text 
TableData = getTableData();
% If Empty return
if(isempty(TableData))
    return
else
    DataOut.data{1} = TableData;
end
% Get the Expiry Date for later use
DataOut.ExpDates{1} = Get_Exp_Dates();
% Get Expiry Dates that are listed in the website to construct separate
% URLS for each month
NextExpiryURL = Get_Next_Expiry_URL();
if(isempty(NextExpiryURL))
   return
end

% Now read Option Tables of each Expiry month
for ik = 1:length(NextExpiryURL)
    urlText = urlread(NextExpiryURL{ik});
    DataOut.ExpDates{ik+1} = Get_Exp_Dates();
    DataOut.data{ik+1} = getTableData();
end
% Clean Up
% Convert the strings into numbers 
f = @(x)[x(:,1) num2cell(str2double(x(:,[2:8]))) x(:,9) num2cell(str2double(x(:,10:end)))];
DataOut.data = cellfun(f,DataOut.data,'uni',false);

goodDataIdx = (~cellfun('isempty',DataOut.data));
DataOut.data = DataOut.data(goodDataIdx );
DataOut.ExpDates = DataOut.ExpDates(goodDataIdx );
% Separate the data into Calls, Puts, Headers
DataOut.Calls = cellfun(@(x) x(:,[1 8 2:7]),DataOut.data,'uni',false);
DataOut.Puts = cellfun(@(x) x(:,[9 8 10:end]),DataOut.data,'uni',false);
DataOut.CPHeaders = {'Symbol','Strike','Last','Change','Bid','Ask','Volume','Open Int'};
DataOut.Headers = {'Symbol','Last','Change','Bid','Ask','Volume','Open Int','Strike',...
    'Symbol','Last','Change','Bid','Ask','Volume','Open Int'};
DataOut.FullOptionData = [DataOut.Headers ; cat(1,DataOut.data{:})];
% Get the Last Price
DataOut.Last = str2num(urlread(['http://download.finance.yahoo.com/d/quotes.csv?s=' symbolid '&f=l1&e=.csv']));

%% Get_Next_Expiry_URL
    function NextExpiry = Get_Next_Expiry_URL()
        % Get the start and end indices and look for a particular text
        Start = regexp(urlText,'View By Expiration:','end');
        end1 = regexp(urlText,'Return to Stacked View...','start');
        
        Data = urlText(Start:end1);
        Data=Data(2:end);
        % Trim the data
        Data=strtrim(Data);
        % Split the data into new lines
        newlines = regexp(Data, '[^\n]*', 'match');
        expr = '<(\w+).*?>.*?</\1>';
        if(isempty(newlines))
            NextExpiry = {};
            return
        end
        % Get the matches of particular expression
        [tok mat] = regexp(newlines{1}, expr, 'tokens', 'match');
        id1= regexp(mat{1},'</b>','start')-1;
        month{1} = mat{1}(4:id1);
        %Month and Next Expiries
        for j = 2:length(mat)-1
            id2 = regexp(mat{j},'">','end');
            id3 = regexp(mat{j},'</a','start');
            if(isempty(id3))
                return
            end
            month{j} = mat{j}(id2+1:id3-1);
            id4 = regexp(mat{j},'"','start');
            NextExpiry{j-1} = ['http://finance.yahoo.com' mat{j}(id4(1)+1:id4(2)-1)]; %#ok<*AGROW>
            NextExpiry{j-1} = regexprep(NextExpiry{j-1},'amp;','');
        end
        
    end
%% Get_Exp_Dates

    function ExpDates = Get_Exp_Dates()
        
        id1 = regexp(urlText,'Options Expiring','end');
        id2 = regexp((urlText(id1+1:id1+51)),'</b>','start');
        ExpDates = strtrim(urlText(id1+1:id1+1+id2-2));
        ExpDates=datestr(datenum(ExpDates,'dddd, mmmm dd,yyyy'));
    end

%% getTableData
    function out = getTableData()
        Main_Pattern = '.*?</table><table[^>]*>(.*?)</table';
        Tables = regexp(urlText, Main_Pattern, 'tokens');
        out = {};
        if(isempty(Tables))
            return
        end
        try
        for TableIdx = 1 : length(Tables)
            
            %Establish a row index
            rowind = 0;
            
            
            % Build cell aray of table data
            
                rows = regexpi(Tables{TableIdx}{:}, '<tr.*?>(.*?)</tr>', 'tokens');
                for rowsIdx = 1:numel(rows)
                    colind = 0;
                    if (isempty(regexprep(rows{rowsIdx}{1}, '<.*?>', '')))
                        continue
                    else
                        rowind = rowind + 1;
                    end
                    
                    headers = regexpi(rows{rowsIdx}{1}, '<th.*?>(.*?)</th>', 'tokens');
                    if ~isempty(headers)
                        for headersIdx = 1:numel(headers)
                            colind = colind + 1;
                            data = regexprep(headers{headersIdx}{1}, '<.*?>', '');
                            if (~strcmpi(data,'&nbsp;'))
                                out{rowind,colind} = strtrim(data);
                            end
                        end
                        continue
                    end
                    cols = regexpi(rows{rowsIdx}{1}, '<td.*?>(.*?)</td>', 'tokens');
                    for colsIdx = 1:numel(cols)
                        if(rowind==1)
                            if(isempty(cols{colsIdx}{1}))
                                continue
                            else
                                colind = colind + 1;
                            end
                        else
                            colind = colsIdx;
                        end
                        % The following code is required to get the sign
                        % of the change in Bid ask prices
                        data = regexprep(cols{colsIdx}{1}, '&nbsp;', ' ');
                        down=false;
                        % If Down is found then it is negative
                        if(~isempty(regexp(data,'"Down"', 'once')))
                            down=true;
                        end
                        data = regexprep(data, '<.*?>', '');
                        if(down)
                            data = ['-' strtrim(data)];
                        end
                        if (~isempty(data))
                            out{rowind,colind} = strtrim(data) ;
                        end
                    end % colsIdx
                end
                
                
        end
        out = out(3:end,:);
        catch %M  %#ok<CTCH> This depends on which version of matlab you are using
               %M.stack
        end
    end
end

MATLAB Built-in Financial Functions

I wrote a Black-Scholes option price calculator function in MATLAB and was working on an implied volatility calculator when I came across some built-in MATLAB functions such as blsprice (Black-Scholes price) and blsimpv (to calculate implied volatility using the Black-Scholes equation).

Pretty cool. I'm sure there are countless others as well...

Put-Call Parity Calculator

If you enjoy the following, consider signing up for the Gödel's Market Newsletter.

I've been reading "Numerical Methods in Finance and Economics, A MATLAB Based Introduction 2nd Edition". They obviously go into put-call parity early on. It's an extremely fascinating concept that given the risk free interest rate (r), the distance until expiration (t), the initial stock price (s), and the strike price (k), that you can then solve for the price of the call (c) given you have the price of the put (p) or vice versa.

No arbitrage financial pricing is extremely interesting (and practical). Here's a MATLAB script I wrote that you can use to look for "arbitrage" opportunity in the put-call parity sense. It'll ask you for several inputs (r,t,s,k,c,p) and tell you what the "call side value" (c + k*e(-r*t)) is and what the "put side value" (p + s) is. You can then compare.

(Example with current SPY options; first part is entry, second is output:)

 % arbitrage; put-call parity checker  
   
 %risk free interest rate  
 r = input('What is the risk free rate? ');  
   
 %time period  
 t = input('How long until expiration? ');  
   
 %stock price, initially  
 s = input('What is the initial stock price? ');  
   
 %strike price  
 k = input('What is the strike price? ');  
   
 %grab call price  
 c = input('What is the call price? ');  
   
 %grab put price  
 p = input('What is the put price? ');  
   
 fprintf('Risk free rate: %d \n', r);  
 fprintf('Time period: %d \n', t);  
 fprintf('Initial stock value: %d \n', s);  
 fprintf('Strike price: %d \n', k);  
 fprintf('Call price: %d \n', c);  
 fprintf('Put price: %d \n', p);  
   
 call_side_value = c + k*exp(-r*t)  
 put_side_value = p + s  

(Edit: fixed an error in call_side_value calculation. Should be good to go.)

(If you've enjoyed this article, consider signing up for the Gödel's Market Newsletter.)

Matlab Stock Correlation Matrix

If you enjoy the following, consider signing up for the Gödel's Market Newsletter.

So, maybe you used the Python Yahoo! data importer script I posted; or, maybe, you have your own price data. And, perhaps you want to run a script to see what the correlation is among the stocks in your database. This code will do just that.

You need to download the MySQL JDBC connector and go through the hassle of adding the .jar file to your Matlab java path. It's worth the struggle, though.

Here's the code with some comments to hopefully help you along (and remind me what I did later):

   
 %test variables to limit matrix sizes...  
 num_stocks = 500;  
 num_close_values = 20;  
 num_log_values = 19;  
   
 %connect to the database  
 conn = database('stocks','root','*****','com.mysql.jdbc.Driver','jdbc:mysql://localhost:3306/stocks');  
   
 %query the symbols in the database  
 str_query_symbols = 'SELECT distinct symbol FROM stock_prices_day';  
 symbol_cell = fetch(conn, str_query_symbols);  
   
 %create matrices to prevent matrix resizing when inserting cell values..  
 mat_log_change = zeros(num_log_values,num_stocks);  
 mat_values_adj_close = zeros(num_close_values,num_stocks);  
   
 %grab adj_close info for a each symbol; then cell2mat it into a larger  
 %matrix  
 for i = 1:num_stocks  
   symbol = cell2mat(symbol_cell(i));  
   str_query_adj_close = ['SELECT adj_close FROM stock_prices_day WHERE symbol = ''' symbol ''' and date > ''2012-06-01'''];  
   mat_values_adj_close(:,i) = cell2mat(fetch(conn, str_query_adj_close));  
     
   %use this to calculate natrual log changes as you move along  
   for j = 1: num_log_values  
     %mat_values_adj_close(j,i)  
     mat_log_change(j,i) = log((mat_values_adj_close(j,i))/(mat_values_adj_close((j+1),i)));  
   end  
 end  
   
 %next calculate correlation coefficients among the various columns  
 mat_correlations = zeros(num_stocks,num_stocks);  
 for i = 1:num_stocks  
   for j = 1:num_stocks  
     A = mat_log_change(:,i);  
     B = mat_log_change(:,j);  
     R = corrcoef(A,B);  
     mat_correlations(i,j) = R(2);  
   end  
 end  
   
 mat_correlations  

(If you've enjoyed this article, consider signing up for the Gödel's Market Newsletter.)