This package calculates the European put and call option prices using the Corrado and Su (1996) model. This method explicitly allows for excess skewness and kurtosis in an expanded Black-Scholes option pricing formula. The approach adapts a Gram-Charlier series expansions of the standard normal density function to yield an option price formula that is the sum of a Black–Scholes option price plus adjustment terms for nonnormal skewness and kurtosis (Corrado and Su, 1997).

For skewness = 0 and kurtosis = 3, the Corrado-Su option prices are equal to the prices obtained using the Black and Scholes (1973) model.

You can download the Matlab code at Corrado and Su (1996) European Option Prices.

References:

Corrado, C.J., and Su T. Skewness and kurtosis in S&P 500 Index returns implied by option prices. Financial Research 19:175–92, 1996.

Corrado, C.J., and Su T. Implied volatility skews and stock return skewness and kurtosis implied by stock option prices. European Journal of Finance 3:73–85, 1997.

Hull, J.C., "Options, Futures, and Other Derivatives", Prentice Hall, 5th edition, 2003.

Luenberger, D.G., "Investment Science", Oxford Press, 1998.

Tags - black scholes , skewness , kurtosis , option

Regime Switching Model with constant transition probability matrix.

Click here for an introduction paper and Matlab codes are here.

Tags - markov , matlab , regime

A key ingredient in these research areas is proper and clean (historical and up-to-date) intra-daily data. On the web there are various resources available, but most of them require a relatively high fee. Other solutions require the use of a specific software. However, there are ways to retrieve intra-daily data for free using Google Finance and also without any software.

If you are familiar with MatLab you can use parts of the package 'Volume Weighted Average Price from Intra-Daily Data' by Semin Ibisevic referenced at Qoppa Investment Society. This package allows you to

(1) retrieve intra-daily stock price data from Google Finance; (2) calculate the VWAP at the end of each trading day; and (3) transform intra-daily data to a daily format. It is a relatively flexible function as it only requires the user to input the ticker symbol and the exchange where the security is listed on. Additionally, the user can define the frequency of the data (1 second or higher) and the period (for instance past 10 days).

If you don't have Matlab you can replicate parts of the code manually. The package above connects with Google Finance and downloads a spreadsheet of intra-daily data consisting of the prices high, open, low, close and volume from: Google Finance get price 1.

You can adjust this to your own preferences by 'seeing' the address as: Google Finance get price 2, where

Hint: to track these inputs, for instance for the Dow Jones Industrial Average, you search the security of interest at Google Finance and then you can find at the top: (INDEXDJX:.DJI) which obviously refers to (EXCHANGE:TICKER).

In a world with increasing competition, such solutions are an handy tool to easier the life of a researcher.

Scholtus, Martin L. and Van Dijk, Dick J. C., High-Frequency Technical Trading: The Importance of Speed (February 28, 2012). Tinbergen Institute Discussion Paper 12-018/4. Available at SSRN: http://papers.ssrn.com/sol3/papers.cfm?abstract_id=2013789.

Ibisevic, Semin. Volume Weighted Average Price from Intra-Daily Data, 2012, MatLab File Exchange. Available online at: http://www.mathworks.com/matlabcentral/fileexchange/36115-volume-weighted-average-price-from-intra-daily-data.

Tags - data , historical , google , intraday

M files:

more can be downloaded at here.

Tags - simulation , sde

Features of the package:

- Support for univariate and multivariate models.

- Support of any number of states and any number of explanatory variables.

- Estimation, by maximum likelihood, of any type of switching setup for the model. This means that you can choose which coefficients in the model, including distribution parameters, are switching states over time.

- A wrapper function for the estimation of regime switching autoregressive models, including multivariate case (MS-VAR) is included in the package.

- The values of standard error for the estimated coefficients can be calculated with 2 different methods.

- Includes a C version of hamilton’s filter that may be used for speeding up the estimation function (see pdf for details).

- Possibility of three distinct distribution assumptions for residual vector (Normal, t or GED).

- Support for reduced/constrained estimation (see pdf document for details).

For instance, as demonstrated here, 3 regimes are simulated as bull, bear and bull markets as

bull1 = normrnd( 0.10, 0.15, 100, 1);

bear = normrnd(-0.01, 0.20, 100, 1);

bull2 = normrnd( 0.10, 0.15, 100, 1);

returns = [bull1; bear; bull2];

bear = normrnd(-0.01, 0.20, 100, 1);

bull2 = normrnd( 0.10, 0.15, 100, 1);

returns = [bull1; bear; bull2];

it is not easy to tell from the return series graph below whether there is regime switch or if yes, when or which part is bear

it is, however, possible by looking at the output of the package, (again, here is one of the reasons I don't buy it, I have to specify how many regimes we expect, how can I know it?)

the graph shows clearly when the regime switching starts.

Overall the package is excellent if you are a lover of the model, check it out yourself. For GAUSS users, I shared a regime switching model library in Gaussbefore.

Tags - regime , switch

For instance, data = get_yahoo_quote({'MSFT','IBM','GOOG','GE'}) returns

1x4 struct array with fields:

symbol

desc

lastTrade

lastTradeTime

dividendDate

open

lastClose

symbol

desc

lastTrade

lastTradeTime

dividendDate

open

lastClose

Tags - data

Nice.

Tags - blog

Therefore I re-introduce a Matlab-only file for

Tags - correlation

Visit Quantitativefinance.co.uk for detail.

Tags - quant

Thanks for all your reply and attention, testing process is fun.

Tags - nag , matlab , fsolve

To make good use of fsolve, the example is to solve the Markowitz problem by finding an optimal portfolio with minimum variance for a targeted return, mathematically, for a portfolio of n risky assets we want to find the solution to:

subject to

here r-bar is a fixed pre-desired level for expected rate of return, and a solution is any portfolio that minimizes the objective function (variance) and offers expected rate r-bar. This is an example of what is called a quadratic program, an optimization problem with a quadratic objective function, and linear constraints. By introducing Lagrange multipliers and taking the first derivative a solution to our Markowitz problem is found by finding a solution to the set of n + 2 linear equations

In the end, the problem falls into the standard framework of linear algebra, and amounts to computing the inverse of a matrix: solve Ax = b, and finally fsolve can be used directly.

To make the problem bigger, I set up a pool with 500 assets, randomly simulate 5-year return and optimize my portfolio under a given target return. 500 assets may be too many for most applications but the set-up of this example is fair enough as we usually have to impose other constraints in reality such as the maximum or minimum weights allocated to a certain industry, etc.

function MarkowitzWeight = fsolve_markowitz_MATLAB()

RandStream.setDefaultStream (RandStream('mt19937ar','seed',1));

% simulate 500 stocks return

SimR = randn(1000,500);

r = mean(SimR); % use mean return as expected return

targetR = 0.02;

rCOV = cov(SimR); % covariance matrix

NumAsset = length(r);

options=optimset('Display','off');

startX = [1/NumAsset*ones(1, NumAsset), 1, 1];

x = fsolve(@fsolve_markowitz,startX,options, r, targetR, rCOV);

MarkowitzWeight = x(1:NumAsset);

end

function F = fsolve_markowitz(x, r, targetR, rCOV)

NumAsset = length(r);

F = zeros(1,NumAsset+2);

weight = x(1:NumAsset); % for asset weights

lambda = x(NumAsset+1);

mu = x(NumAsset+2);

for i = 1:NumAsset

F(i) = sum(weight.*rCOV(i,:))-lambda*r(i)-mu;

end

F(NumAsset+1) = sum(weight.*r)-targetR;

F(NumAsset+2) = sum(weight)-1;

end

RandStream.setDefaultStream (RandStream('mt19937ar','seed',1));

% simulate 500 stocks return

SimR = randn(1000,500);

r = mean(SimR); % use mean return as expected return

targetR = 0.02;

rCOV = cov(SimR); % covariance matrix

NumAsset = length(r);

options=optimset('Display','off');

startX = [1/NumAsset*ones(1, NumAsset), 1, 1];

x = fsolve(@fsolve_markowitz,startX,options, r, targetR, rCOV);

MarkowitzWeight = x(1:NumAsset);

end

function F = fsolve_markowitz(x, r, targetR, rCOV)

NumAsset = length(r);

F = zeros(1,NumAsset+2);

weight = x(1:NumAsset); % for asset weights

lambda = x(NumAsset+1);

mu = x(NumAsset+2);

for i = 1:NumAsset

F(i) = sum(weight.*rCOV(i,:))-lambda*r(i)-mu;

end

F(NumAsset+1) = sum(weight.*r)-targetR;

F(NumAsset+2) = sum(weight)-1;

end

I performed the calculations tic; fsolve_markowitz_MATLAB(); toc on a computer running Windows XP, MATLAB 2008b and got the following timings (averaged over 10 runs): 3.35 seconds.

So Michael, are you able to speed up my example for 19.7 times faster?

Tags - nag , matlab , markowitz , fsolve

You can view the agenda and register online at http://www.mathworks.com/company/events/conferences/mc2010-linkedin/index.html, a few possible interesting events are:

This master class will show how you can take advantage of symbolic calculations within the MATLAB environment for modelling, simulation, and design tasks. Using the latest features in Symbolic Math Toolbox™, including the new MuPAD notebook interface, two case studies will be examined. A wind turbine model will be developed, documented, and integrated with MATLAB as part of a design optimization study. We will also examine the use of MathWorks tools for symbolic and numeric calculations using an example from the BLOODHOUND Project, which involves a World Land Speed Record attempt and aims to inspire young people to pursue careers in science, technology, engineering, and mathematics.

In this master class, you'll see how you can use MATLAB and MathWorks parallel computing tools to solve computationally and data-intensive problems taking advantage of recent advances in computer hardware, from multiprocessor machines to computer clusters. You will learn how to utilise multiple cores in your desktop machine through the new capabilities of MATLAB and Parallel Computing Toolbox. We will also show how to scale applications to computer clusters without changing the code.

FTSE Group is a world leader in the creation and management of over 120,000 equity, bond, and alternative asset class indices. Index creation is FTSE's sole business, and the demand for new indices to support more sophisticated investment strategies is growing. FTSE uses MATLAB to enable analysts to work closely with clients while developing algorithms for customised index solutions.

This presentation will describe the results of a proof-of-concept project to demonstrate how MATLAB can be used in a grid to maintain and improve service delivery times while satisfying demand and shortening product development life cycles.

Tags - matlab , conference

Quotation

Nonlinear principal component analysis (NLPCA) is commonly seen as a nonlinear generalization of standard principal component analysis (PCA). It generalizes the principal components from straight lines to curves (nonlinear). Thus, the subspace in the original data space which is described by all nonlinear components is also curved.

Nonlinear PCA can be achieved by using a neural network with an autoassociative architecture also known as autoencoder, replicator network, bottleneck or sandglass type network. Such autoassociative neural network is a multi-layer perceptron that performs an identity mapping, meaning that the output of the network is required to be identical to the input. However, in the middle of the network is a layer that works as a bottleneck in which a reduction of the dimension of the data is enforced. This bottleneck-layer provides the desired component values (scores).

Nonlinear PCA can be achieved by using a neural network with an autoassociative architecture also known as autoencoder, replicator network, bottleneck or sandglass type network. Such autoassociative neural network is a multi-layer perceptron that performs an identity mapping, meaning that the output of the network is required to be identical to the input. However, in the middle of the network is a layer that works as a bottleneck in which a reduction of the dimension of the data is enforced. This bottleneck-layer provides the desired component values (scores).

Interested readers shall download the

Tags - pca , toolbox

Let's say you have a csv file with structure like below

How to import it in Matlab? csvread() obviously doesn't work as this file has mixed format, and it returns error

Quotation

??? Error using ==> dlmread at 145

Mismatch between file and format string.

Trouble reading number from file (row 1, field 1) ==> CUSIP

Error in ==> csvread at 52

m=dlmread(filename, ',', r, c);

Mismatch between file and format string.

Trouble reading number from file (row 1, field 1) ==> CUSIP

Error in ==> csvread at 52

m=dlmread(filename, ',', r, c);

Luckily I found a function called

Tags - matlab , csv

Command-line access to MATLAB running on your computer

Access to your MATLAB workspace

Ability to view MATLAB figures on your iPhone

Record of commands typed on the iPhone in your command history

Custom keyboard

Ability to connect to MATLAB running on Windows, Mac, and Linux

Download the free app

Tags - matlab , iphone , gphone

Send Text Message to Cell Phone is such a great file I found recently, bascially what it does is to send email via sendmail function of Matlab from your gmail box to your cell phone carrier, and then your cell phone carrier forwards the email to you as a text message.

However, it works for US based cell phones only, I have tried on my UK T-mobile phone and it seems UK T-mobile doesn't support such a mail to SMS service (correct me if I am wrong).

Fortunately, I came across SMS service website which allows people to send up to 3 free email to SMS per day, it should be enough for our use in Matlab. Add the following line in the switch case after line 55

case 'uk'; emailto = strcat(number,'@x-onsms.com');

that's it, the email will be delivered as an SMS to your mobile. Do let me know if you are aware of a better alternative.

So what you need to do is to put the function send_text_message at the end of your file, it will then send you a message automatically, for example

send_text_message('079-123-456','UK', 'Desk 12 Calculation Done','Now you can shut down the computer')

What else can it be used? stock price alert? profit threshold alarm?

Possible error

Depends on your Matlab version and firewall setting, you may notice the following errors:

530 5.7.0 Must issue a STARTTLS command first

This could happen for MATLAB 7.1 (R14SP3) and before, you may have to upgrade your version.

Connection timed out: connect

This is due to your firewall or anti-virus software setting. you are not allowed to send email from port 25. What you shall do is too add an exception and let your computer know this action is safe. For instance, I have McAfee, to add an exception, open its

add Matlab.exe as a process to be excluded, save it, done.

Ros, you don't have to check computers one by one. Sounds useful? download the file at http://www.mathworks.com/matlabcentral/fileexchange/16649, don't forget to change the email address and password at the beginning of the file.

Tags - matlab , sms

Quotation

Writing program code is a good way of debugging your thinking - Bill Venables

Believe it or not, on average programmers spend 70% of their working time in debugging (I don't know where this number is from but I do believe so). A good debug tool, command or even habit will definitely improve your work efficiency, shorten working hour, and enjoy more the world cup. When it comes to Matlab and R, I have to say the debugging in Matlab is straightforward but is more comprehensive in R. Here is an introductory video demonstrating how to debug and understand Matlab codes:

NEWFCN is a good function I always use when creating new files, basically it creates a M-File having the entered filename and a specific structure which helps for creating the main function structure. The actual working MATLAB Version will be also captured. For example, running

function testnewfcn()

% TESTNEWFCN ...

%

% ...

%% AUTHOR : aBiao @ mathfinance.cn

%% $DATE : 28-May-2010 18:55:27 $

%% $Revision : 1.00 $

%% DEVELOPED : 7.1.0.246 (R14) Service Pack 3

%% FILENAME : testnewfcn.m

disp(' !!! You must enter code into this file < testnewfcn.m > !!!')

% ===== EOF ====== [testnewfcn.m] ======

% TESTNEWFCN ...

%

% ...

%% AUTHOR : aBiao @ mathfinance.cn

%% $DATE : 28-May-2010 18:55:27 $

%% $Revision : 1.00 $

%% DEVELOPED : 7.1.0.246 (R14) Service Pack 3

%% FILENAME : testnewfcn.m

disp(' !!! You must enter code into this file < testnewfcn.m > !!!')

% ===== EOF ====== [testnewfcn.m] ======

Did I mention you can customize the exact style by modifying newfcn.m to the appearance you like? just that simple!

Download at http://www.mathworks.com/matlabcentral/fileexchange/6408

Tags - matlab

Quotation

FinMetrics is MATLAB based, open source quantitative portfolio management environment. Built on concepts of bottom-up approach to application design, it allows users to define most basic, low level building blocks, e.g. assets and transactions, to be further pieced together in a higher level objects, e.g. positions or portfolios. Data analysis and statistics function, implemented within the environment and native to MATLAB, enable users to conduct scenario analysis, stress testing, performance measurement and attribution, risk measurement and attribution, design hedge strategies, etc. Open architecture of the environment allows users to work with objects of any level, depending on their requirements and expertise. The object structure and data types are specifically designed to make integration with MATLAB and native FinMetrics functions as easy as possible. FinMetrics user interface application and MATLAB scripting may be utilized to facilitate or automate complex and repetitive tasks, as well as extend functionality of the environment.

Download it at http://www.mathworks.com/matlabcentral/fileexchange/27778-finmetrics if interested.

Tags - matlab

It is all right to do like that, but the whole process becomes extremely simple with the cell mode in Matlab, it generates report automatically for us, any change you make regarding description, codes and results will be updated by clicking a simple button:

The following steps outline the procedure,

First, opening an editor in Matlab by choosing File -> New -> M-file;

Second, selecting cell and enable cell mode in the editor;

Third, inserting cell divider under cell window, it allows you to type a name for each section, similar with the title for chapters in Word;

Fourth, depending on your reports, inserting Text Markup, where you are able to insert description text, bold text, etc. most importantly, it supports TeX equations;

Finally, clicking "Publish to HTML", a nice-looking report is generated, and if necessary, changing the codes and re-running to update the report.

Second, selecting cell and enable cell mode in the editor;

Third, inserting cell divider under cell window, it allows you to type a name for each section, similar with the title for chapters in Word;

Fourth, depending on your reports, inserting Text Markup, where you are able to insert description text, bold text, etc. most importantly, it supports TeX equations;

Finally, clicking "Publish to HTML", a nice-looking report is generated, and if necessary, changing the codes and re-running to update the report.

For a simple demonstration, I write the files as follows.

%% Cell mode publish demonstration

% abiao @ mathfinance.cn

%% Using Black Scholes formula as an example

% The value of a call and put options in terms of the Black–Scholes parameters

% is:

%

% $$C(S,t) = SN(d_1) - Ke^{-r(T - t)}N(d_2) \,$$

%

% $$P(S,t) = Ke^{-r(T-t)} - S + (SN(d_1) - Ke^{-r(T - t)}N(d_2)) = Ke^{-r(T-t)} - S + C(S,t)$$

%

% $$d_1 = \frac{\ln(\frac{S}{K}) + (r + \frac{\sigma^2}{2})(T - t)}{\sigma\sqrt{T - t}}$$

%

% $$d_2 = d_1 - \sigma\sqrt{T - t}$$

%% Sample code, for simplicity, I use the embedded command

[call, put] = blsprice(10,9,0.02,2,0.3);

fprintf('Call option value is %g. \n',call)

fprintf('Put option value is %g. \n',put)

call1 = [];

put1 = [];

for i = 1:10

call1(i) = blsprice(5+i,9,0.02,2,0.3);

end

%% Sample plot for demonstration only

plot(5+(1:10), call1)

% abiao @ mathfinance.cn

%% Using Black Scholes formula as an example

% The value of a call and put options in terms of the Black–Scholes parameters

% is:

%

% $$C(S,t) = SN(d_1) - Ke^{-r(T - t)}N(d_2) \,$$

%

% $$P(S,t) = Ke^{-r(T-t)} - S + (SN(d_1) - Ke^{-r(T - t)}N(d_2)) = Ke^{-r(T-t)} - S + C(S,t)$$

%

% $$d_1 = \frac{\ln(\frac{S}{K}) + (r + \frac{\sigma^2}{2})(T - t)}{\sigma\sqrt{T - t}}$$

%

% $$d_2 = d_1 - \sigma\sqrt{T - t}$$

%% Sample code, for simplicity, I use the embedded command

[call, put] = blsprice(10,9,0.02,2,0.3);

fprintf('Call option value is %g. \n',call)

fprintf('Put option value is %g. \n',put)

call1 = [];

put1 = [];

for i = 1:10

call1(i) = blsprice(5+i,9,0.02,2,0.3);

end

%% Sample plot for demonstration only

plot(5+(1:10), call1)

The final report is shown at http://www.mathfinance.cn/attachment/CellPublishDemo.html, looks fantastic.

Tags - matlab

Needless to say, comment is crucial for programming, it helps ourselves to debug our thoughts and other colleagues to understand the codes efficiently, which is especially indispensable as each programmer has his own coding style. However, under certain circumstances

% test the comment strip command using mlstripcommentsfile(infile, outfile)

% by abiao @ mathfinance.cn

a = rnorm(10,1); % check if you can remove me

% one more line

b = rnorm(5,5);

c = rnorm(3,3)...

+ rnorm(3,3); %now i am here

% test finished

% by abiao @ mathfinance.cn

a = rnorm(10,1); % check if you can remove me

% one more line

b = rnorm(5,5);

c = rnorm(3,3)...

+ rnorm(3,3); %now i am here

% test finished

after running the command mlstripcommentsfile('test.m', 'test-client.m'), I got the file test-client.m to send to others

a = rnorm(10,1);

b = rnorm(5,5);

c = rnorm(3,3)...

+ rnorm(3,3);

Nice, isn't it? you can download the toolbox at http://www.mathworks.com/matlabcentral/fileexchange/4645.

PS: the original codes use old Matlab version so you will get the warning "Warning: The 'tokenize' option for regexprep is now the default behavior." You may choose to modify the codes by deleting 'tokenize'.

Tags - matlab

I skip the technical section, which can be found in the original paper Beyond Black-Litterman in Practice: A Five-Step Recipe to Input Views on Non-Normal Markets, and the Matlab codes can be downloaded at http://www.mathworks.com/matlabcentral/fileexchange/9061, in the Extra->COP folder. Suppose we would like to invest in the US treasury market at a weekly investment horizon, and we are interested into the following six key interest rates: 6 month, 2 year, 5 year, 10 year, 20 year and 30 year. For illustration, we use Monte Carlo simulation to generate 100,000 scenarios based on a t Copula with skew t marginal distribution, the sample mean, standard deviation, skewness and kurtosis of the six key rates are shown in table 1.

(table 1)

We can easily see from the table 1 that all of the key rates, in particular the short period rates, are non-normally distributed as the kurtosis is significantly larger than 3, which corresponds to normal distribution.

Compared with table 1, we notice the means of 2y and 10y are decreased significantly, and the means of 5y and 20y are increased, which are consistent with our assumed views, since we are bearish on 2y and 10y, and bullish on 5y and 20y.

as expected, the means of 2y and 10y are increased significantly, the means of 5y and 20y are increased, and the means of 6m and 30y are again the same as those of prior market.

identical with what we expected, when compared with the results of case 1, the means of 2y and 10y are decreased more significantly, and the means of 5y and 20y are increased more pronounced. For all views, the even moments are only marginally affected. Our robust analysis demonstrates the strong consistence of the COP approach.

Based on the posterior market values we are able to do return mapping and portfolio optimization.

These two posts are a short summary of the original paper and Matlab implementation, please read the source for detail.

Tags - black-litterman

Markowitz mean-variance model is the foundation of modern portfolio theory, it tries to maximize return and minimize risk by carefully choosing different assets. Although Markowitz model is widely used in practice in the financial industry, it has at least the following shortcomings. First, in principle the Markowitz model offers a solution once the expected returns and covariance of the assets are known. Although the covariance of a few assets can be adequately estimated, it is difficult to come up with reasonable estimates of expected returns. Initial estimates often lead to extreme portfolio solutions and have to be adjusted many times before they can be considered acceptable by the decision makers. Second, Markowitz model assumes the future distribution of assets returns will be exactly the same as that of historical returns, therefore it does not allow an investor to have views on the direction of assets growth.

Black-Litterman model bypasses these problems by not requiring the user to input estimates of expected return; instead it assumes that the initial expected returns are whatever is required so that the equilibrium asset allocation is equal to what we observe in the markets. The user is only required to state how his assumptions about expected returns differ from the market's and to state his degree of confidence in the alternative assumptions. From this, the Black-Litterman model computes the desired mean-variance efficient asset allocation. It is right because of this character that makes Black-Litterman model be popular both in industry and academia since it was initially developed in 1990.

However, the original Black-Litterman model helps portfolio managers to compute the distribution of the posterior market that incorporates their subjective views with respect to the prior distribution. All these assumptions are being made considering the entire scenario is normally distributed. In fact, it is frequently observed that returns in equity and other markets are not normally distributed. At the same time, you would expect a practitioner to input his views in a less informative manner. For instance, instead of a view like ‘I am 90% sure that stock A will have a 2% expected return’, he may have a more accurate view such as ‘I am 90% sure that the return of stock A will stay between -1% and 1% with equal chance’, or for an Asian option Trader, his views may be a set of prices at several specific monitoring times, when the option is supposed to be exercised.

Meucci (2006) solves the above-mentioned issues of the original Black-Litterman model by developing an approach called Copula-Opinion Pooling (COP), the author relaxes the normal distribution assumption and in principle allows any distribution. However, the proposed implementation in the paper has a distinct distribution under certain assumptions. Indeed, the market has been assumed to be skew-t distributed, this strong assumption makes the application of the approach less evident, although the implementation is pretty straightforward for any distribution by Monte Carlo simulation, because Monte Carlo simulation can be applied for almost any distribution and is suitable for any dimension (at least in theory), COP approach can handle a portfolio with any number of risk factors.

The second part will be about the Matlab implementation and results analysis.

Tags - black-litterman , allocation

those different formats really bring a problem, as I need to match the dates and compare the prices, it is easy to convert a cell format date to what matlab recognizes using cell2mat() and then datenum(), for instance, matlab returns a number 733737 for '25/11/2008', however, how can I compare that with 20081125 then?

what I am thinking is:

It sounds lengthy, do you have a better idea? please share with us by leaving a comment, cheers.

Stanley suggests to use the following code: datenum(datevec(num2str(dates),'yyyymmdd')), here dates is double number, perfect, it is much better than my method as it allows us not only for date comparison but also for calculation, what's more, it is more efficient. many thanks, stanley.

Tags - matlab , excel

Specifically, Kalman Filter is applied to estimate the parameters of a Cox Ingersoll Ross (CIR) one factor interest rate model, (Vasicek model is simplier than CIR, so the latter is chosen as an example), it is a widely used mean-reverting process with SDE

A three-factor CIR model has a measurement equation

and a transition equation

source from the paper "affine term structure models: theory and implementation" downloaded at www.bankofcanada.ca/en/res/wp/2001/wp01-15a.pdf

I skip the derivation part and recommend the following two papers: "estimating and testing exponential-affine term structure models by kalman filter " and "affine term structure models: theory and implementation" to understand the transition and measurement equations. Below are the sample Matlab implementation:

function [para, sumll] = TreasuryYieldKF()

% author: biao from www.mathfinance.cn

%%CIR parameter estimation using Kalman Filter for given treasury bonds yields

% check paper ""estimating and testing exponential-affine term structure

% models by kalman filter " and "affine term structure models: theory and

% implementation" for detail

% S(t+1) = mu + F S(t) + noise(Q)

% Y(t) = A + H S(t) + noise(R)

% read data Y

Y = xlsread('ir.xls');

[nrow, ncol] = size(Y);

tau = [1/4 1/2 1 5]; % stand for 3M, 6M, 1Y, 5Y yield

para0 = [0.05, 0.1, 0.1, -0.1, 0.01*rand(1,ncol).*ones(1,ncol)];

[x, fval] = fmincon(@loglik, para0,[],[],[],[],[0.0001,0.0001,0.0001, -1, 0.00001*ones(1,ncol)],[ones(1,length(para0))],[],[],Y, tau, nrow, ncol);

para = x;

sumll = fval;

end

% author: biao from www.mathfinance.cn

%%CIR parameter estimation using Kalman Filter for given treasury bonds yields

% check paper ""estimating and testing exponential-affine term structure

% models by kalman filter " and "affine term structure models: theory and

% implementation" for detail

% S(t+1) = mu + F S(t) + noise(Q)

% Y(t) = A + H S(t) + noise(R)

% read data Y

Y = xlsread('ir.xls');

[nrow, ncol] = size(Y);

tau = [1/4 1/2 1 5]; % stand for 3M, 6M, 1Y, 5Y yield

para0 = [0.05, 0.1, 0.1, -0.1, 0.01*rand(1,ncol).*ones(1,ncol)];

[x, fval] = fmincon(@loglik, para0,[],[],[],[],[0.0001,0.0001,0.0001, -1, 0.00001*ones(1,ncol)],[ones(1,length(para0))],[],[],Y, tau, nrow, ncol);

para = x;

sumll = fval;

end

flip over to next page...

function sumll = loglik(para,Y, tau, nrow, ncol) %calculate log likelihood

% initialize the parameter for CIR model

theta = para(1); kappa = para(2); sigma = para(3); lambda = para(4);

%volatility of measurement error

sigmai = para(5:end);

R = eye(ncol);

for i = 1:ncol

R(i,i) = sigmai(i)^2;

end

dt = 1/12; %monthly data

initx = theta;

initV = sigma^2*theta/(2*kappa);

% parameter setting for transition equation

mu = theta*(1-exp(-kappa*dt));

F = exp(-kappa*dt);

% parameter setting for measurement equation

A = zeros(1, ncol);

H = A;

for i = 1:ncol

AffineGamma = sqrt((kappa+lambda)^2+2*sigma^2);

AffineBeta = 2*(exp(AffineGamma*tau(i))-1)/((AffineGamma+kappa+lambda)*(exp(AffineGamma*tau(i))-1)+2*AffineGamma);

AffineAlpha = 2*kappa*theta/(sigma^2)*log(2*AffineGamma*exp((AffineGamma+kappa+lambda)*tau(i)/2)/((AffineGamma+kappa+lambda)*...

(exp(AffineGamma*tau(i))-1)+2*AffineGamma));

A(i) = -AffineAlpha/tau(i);

H(i) = AffineBeta/tau(i);

end

%now recursive steps

AdjS = initx;

VarS = initV;

ll = zeros(nrow,1); %log-likelihood

for i = 1:nrow

PredS = mu+F*AdjS; %predict values for S and Y

Q = theta*sigma*sigma*(1-exp(-kappa*dt))^2/(2*kappa)+sigma*sigma/kappa*(exp(-kappa*dt)-exp(-2*kappa*dt))*AdjS;

VarS = F*VarS*F'+Q;

PredY = A+H*PredS;

PredError = Y(i,:)-PredY;

VarY = H'*VarS*H+R;

InvVarY = inv(VarY);

DetY = det(VarY);

%updating

KalmanGain = VarS*H*InvVarY;

AdjS = PredS+KalmanGain*PredError';

VarS = VarS*(1-KalmanGain*H');

ll(i) = -(ncol/2)*log(2*pi)-0.5*log(DetY)-0.5*PredError*InvVarY*PredError';

end

sumll = -sum(ll);

end

% initialize the parameter for CIR model

theta = para(1); kappa = para(2); sigma = para(3); lambda = para(4);

%volatility of measurement error

sigmai = para(5:end);

R = eye(ncol);

for i = 1:ncol

R(i,i) = sigmai(i)^2;

end

dt = 1/12; %monthly data

initx = theta;

initV = sigma^2*theta/(2*kappa);

% parameter setting for transition equation

mu = theta*(1-exp(-kappa*dt));

F = exp(-kappa*dt);

% parameter setting for measurement equation

A = zeros(1, ncol);

H = A;

for i = 1:ncol

AffineGamma = sqrt((kappa+lambda)^2+2*sigma^2);

AffineBeta = 2*(exp(AffineGamma*tau(i))-1)/((AffineGamma+kappa+lambda)*(exp(AffineGamma*tau(i))-1)+2*AffineGamma);

AffineAlpha = 2*kappa*theta/(sigma^2)*log(2*AffineGamma*exp((AffineGamma+kappa+lambda)*tau(i)/2)/((AffineGamma+kappa+lambda)*...

(exp(AffineGamma*tau(i))-1)+2*AffineGamma));

A(i) = -AffineAlpha/tau(i);

H(i) = AffineBeta/tau(i);

end

%now recursive steps

AdjS = initx;

VarS = initV;

ll = zeros(nrow,1); %log-likelihood

for i = 1:nrow

PredS = mu+F*AdjS; %predict values for S and Y

Q = theta*sigma*sigma*(1-exp(-kappa*dt))^2/(2*kappa)+sigma*sigma/kappa*(exp(-kappa*dt)-exp(-2*kappa*dt))*AdjS;

VarS = F*VarS*F'+Q;

PredY = A+H*PredS;

PredError = Y(i,:)-PredY;

VarY = H'*VarS*H+R;

InvVarY = inv(VarY);

DetY = det(VarY);

%updating

KalmanGain = VarS*H*InvVarY;

AdjS = PredS+KalmanGain*PredError';

VarS = VarS*(1-KalmanGain*H');

ll(i) = -(ncol/2)*log(2*pi)-0.5*log(DetY)-0.5*PredError*InvVarY*PredError';

end

sumll = -sum(ll);

end

I also attach the ir.xls file to test the codes, you should get the parameters values as:

theta: 0.0613

kappa: 0.2249

sigma: 0.0700

lambda: -0.1110

Click to download

Tags - filter

In the book

1, known dividend yield. For instance, there will be a 3% dividend 3 months later (3% of the stock price), it is straightforward to handle it as the binomial tree is recombined when the nodes are multiplied by a percentage, so basically what we need to do is to construct a tree like usual before ex-dividend date, and then shift all the left tree nodes down by (1-dividend yield), that's it, the number of nodes are the same as for non-dividend binomial tree;

(source from

2, known dollar dividend. For instance, there will be a 2.5 dollar dividend 3 months later, so before ex-dividend date the binomial tree is constructed as usual but exactly at the date after ex-dividend, the whole nodes are shifted down by 2.5 dollar, and then a new binomial tree is constructed, because the nodes are shifted by an absolute amount number, the new binomial tree is not recombined any more, which means much more nodes than the non-dividend case. Specifically, as pointed by Hull, when i = k+m, there are m(k+2) rather than k+m+1 nodes. The issue becomes more challenging when we increase the number of dividends. Fortunately, there is a simpler way to get around of this difficulty by dividing the stock price into two components: an uncertain part and a part that is the present value of all future dividends during the life of the option. Please check the book for detail Options, Futures, and Other Derivatives, 7th Economy Edition with CD.

(source from

Should you are interested into a sample implementation in Matlab of Binomial Tree Option Pricing with Discrete Dividends, take a look at the file http://www.ualberta.ca/dept/aict/bluejay/usr/local/matlab-6.5/toolbox/finance/finance/binprice.m.

Tags - dividend , option

Found a paper

Take a look if you are interested, in the Appendix the R and Matlab codes are given for a better understanding. http://papers.ssrn.com/sol3/papers.cfm?abstract_id=1566975

Tags - stochastic , heuristic

I happened to find Cleve Moler solved Sudoku using recursive backtracking method and the speed is fast, the minor error is it doesn't report an error when the initial 9*9 matrix you give violates Sudoku rules (you can test it later, the code starts to run and returns a result even you give repetitive numbers along a same row or column). Since building a GUI in matlab isn't easy, I choose to build a

Once you install Excel Link module and turn it on, you will notice the short-cuts on the excel menu bar looking like

intuitively, those buttons stand for "start matlab", "send data to matlab", "retrieve data from matlab", and "execute the matlab command", respectively.

For my case, I first download the Sudoku M code at http://www.mathworks.co.uk/company/newsletters/news_notes/2009/clevescorner.html?s_cid=ACD0210ukTA2&s_v1=8728847_1-BR7DSN, then open an excel file, write down a 9*9 matrix "X"

then define a zone to fetch the calculation result, which should also be 9*9

finally do steps:

<- put data to matlab, MLPutmatrix("X",X)

<- solve the problem in matlab, MLEvalString("X=sudoku(X)")

<- get results from matlab, MLGetMatrix("X","NewX")

result is then retrieved immediately to Sudoku Spreadsheet from matlab

Straightforward to run your Matlab function in excel, isn't it? alternatively you can use Matlab builder for excel.

Click to download

Tags - matlab , excel

1,

2,

3,

4,

5,

6,

Below is sample Matlab codes I wrote for comparision, a single dividend is used for simplicity, results similar to the table listed in the paper

function callprice = DiscreteDividend(s,k,r,t,vol,d,dt)

%compare different methods for a single discrete dividend adjustments,

%read paper by Haug for detail;

%dt: dividend time

BSprice = blsprice(s,k,r,t,vol,0);

AdjS = s-exp(-r*dt)*d;

Escrowed = blsprice(AdjS,k,r,t,vol,0);

%%%%%%%%%Chriss, 1997%%%%%%%%%%

vol1 = vol*s/(s-d*exp(-r*dt));

Chriss = blsprice(AdjS,k,r,t,vol1,0);

%%%%%%%%%Haug, 1998%%%%%%%%%%

vol2 = sqrt(vol^2*dt+vol1^2*(t-dt)/t);

OldHaug = blsprice(AdjS,k,r,t,vol2,0);

%%%%%%%%Bos et al. (2003)%%%%%%%%%%

lns = log(s);

lnk = log((k+d*exp(-r*dt))*exp(-r*t));

z1 = (lns-lnk)/(vol*sqrt(t))+vol*sqrt(t)/2;

z2 = z1+vol*sqrt(t)/2;

vol3 = sqrt(vol^2+vol*sqrt(pi/(2*t))*(4*exp(z1^2/2-lns)*d*exp(-r*dt)*...

(normcdf(z1)-normcdf(z1-vol*dt/sqrt(t)))+exp(z2^2/2-2*lns)*d^2*...

exp(-r*2*dt)*(normcdf(z2)-normcdf(z2-2*vol*dt/sqrt(t)))));

Bos = blsprice(AdjS,k,r,t,vol3,0);

%%%%%%%%%Haug, 2003%%%%%%%%%%%%%%%

NewHaug = exp(-r*dt)*(quad(@(x)blsprice(x-d,k,r,t-dt,vol,0).*lognpdf(x,lns+(r-0.5*vol^2)*dt,vol*sqrt(dt)), d, k+d)...

+quad(@(x)blsprice(x-d,k,r,t-dt,vol,0).*lognpdf(x,lns+(r-0.5*vol^2)*dt,vol*sqrt(dt)), k+d, 20*s));

callprice = [BSprice, Escrowed, Chriss, OldHaug, Bos, NewHaug];

%compare different methods for a single discrete dividend adjustments,

%read paper by Haug for detail;

%dt: dividend time

BSprice = blsprice(s,k,r,t,vol,0);

AdjS = s-exp(-r*dt)*d;

Escrowed = blsprice(AdjS,k,r,t,vol,0);

%%%%%%%%%Chriss, 1997%%%%%%%%%%

vol1 = vol*s/(s-d*exp(-r*dt));

Chriss = blsprice(AdjS,k,r,t,vol1,0);

%%%%%%%%%Haug, 1998%%%%%%%%%%

vol2 = sqrt(vol^2*dt+vol1^2*(t-dt)/t);

OldHaug = blsprice(AdjS,k,r,t,vol2,0);

%%%%%%%%Bos et al. (2003)%%%%%%%%%%

lns = log(s);

lnk = log((k+d*exp(-r*dt))*exp(-r*t));

z1 = (lns-lnk)/(vol*sqrt(t))+vol*sqrt(t)/2;

z2 = z1+vol*sqrt(t)/2;

vol3 = sqrt(vol^2+vol*sqrt(pi/(2*t))*(4*exp(z1^2/2-lns)*d*exp(-r*dt)*...

(normcdf(z1)-normcdf(z1-vol*dt/sqrt(t)))+exp(z2^2/2-2*lns)*d^2*...

exp(-r*2*dt)*(normcdf(z2)-normcdf(z2-2*vol*dt/sqrt(t)))));

Bos = blsprice(AdjS,k,r,t,vol3,0);

%%%%%%%%%Haug, 2003%%%%%%%%%%%%%%%

NewHaug = exp(-r*dt)*(quad(@(x)blsprice(x-d,k,r,t-dt,vol,0).*lognpdf(x,lns+(r-0.5*vol^2)*dt,vol*sqrt(dt)), d, k+d)...

+quad(@(x)blsprice(x-d,k,r,t-dt,vol,0).*lognpdf(x,lns+(r-0.5*vol^2)*dt,vol*sqrt(dt)), k+d, 20*s));

callprice = [BSprice, Escrowed, Chriss, OldHaug, Bos, NewHaug];

For example, the results of a $7 dividend after 0.5 year are (DiscreteDividend(100, 100,0.06,1,0.3,7,0.5)): 14.7171 10.6932 11.5001 11.1039 11.0781 11.1062, respectively.

Reading the original paper

Tags - dividend , option

By building a self-financing portfolio consisting a unit of option at strike K, -delta1 units of underlying asset, and -xi units of options at strike ki, Castagna and Mercurio (2007) calculate the weight xi for three given quotes with the help of Ito's lemma and then approximate the European option value under

Below is a simple Matlab code to price a call option based on Castagna and Mercurio (2007):

% option price under Vanna volga model for any strike k

% sigma2 is ATM implied vol, k2 is ATM strike

s = 1.205;

t = 94/365;

r = -log(0.9902752)/t;

rf = -log(0.9945049)/t;

sigma1 = 0.0979;

sigma3 = 0.0929;

sigma2 = 0.09375;

k1 = 1.172;

k3 = 1.2504;

k2 = 1.2115;

k = 1.24;

Vega1 = Vega(s,k1,r,t,sigma2,rf);

Vega3 = Vega(s,k3,r,t,sigma2,rf);

Vegak = Vega(s,k,r,t,sigma2,rf);

x1 = Vegak*log(k2/k)*log(k3/k)/(Vega1*log(k2/k1)*log(k3/k1));

x3 = Vegak*log(k/k1)*log(k/k2)/(Vega3*log(k3/k1)*log(k3/k2));

price = blsprice(s,k,r,t,sigma2,rf)+x1*(blsprice(s,k1,r,t,sigma1,rf)...

-blsprice(s,k1,r,t,sigma2,rf))+x3*(blsprice(s,k3,r,t,sigma3,rf)...

-blsprice(s,k3,r,t,sigma2,rf));

% sigma2 is ATM implied vol, k2 is ATM strike

s = 1.205;

t = 94/365;

r = -log(0.9902752)/t;

rf = -log(0.9945049)/t;

sigma1 = 0.0979;

sigma3 = 0.0929;

sigma2 = 0.09375;

k1 = 1.172;

k3 = 1.2504;

k2 = 1.2115;

k = 1.24;

Vega1 = Vega(s,k1,r,t,sigma2,rf);

Vega3 = Vega(s,k3,r,t,sigma2,rf);

Vegak = Vega(s,k,r,t,sigma2,rf);

x1 = Vegak*log(k2/k)*log(k3/k)/(Vega1*log(k2/k1)*log(k3/k1));

x3 = Vegak*log(k/k1)*log(k/k2)/(Vega3*log(k3/k1)*log(k3/k2));

price = blsprice(s,k,r,t,sigma2,rf)+x1*(blsprice(s,k1,r,t,sigma1,rf)...

-blsprice(s,k1,r,t,sigma2,rf))+x3*(blsprice(s,k3,r,t,sigma3,rf)...

-blsprice(s,k3,r,t,sigma2,rf));

where Vega is a function to compute vega under black scholes formula

function VegaValue = Vega(s,k,r,t,sigma,rf)

d1 = (log(s/k)+(r-rf+0.5*sigma^2)*t)/(sigma*sqrt(t));

VegaValue = s*exp(-rf*t)*sqrt(t)*normpdf(d1,0,1);

d1 = (log(s/k)+(r-rf+0.5*sigma^2)*t)/(sigma*sqrt(t));

VegaValue = s*exp(-rf*t)*sqrt(t)*normpdf(d1,0,1);

Implied volatilities curve is therefore easily achieved by inverting VV pricing model. Interested ppl please refer to http://www.risk.net/risk/technical-paper/1506580/the-vanna-volga-method-implied-volatilities or an advanced one www.mathfinance.de/wystup/papers/wystup_vannavolga_eqf.pdf

Tags - volatility

Quotation

Abstract: One of the most widely used option valuation procedures among practitioners is a version of Black-Scholes in which implied volatilities are smoothed across strike prices and maturities. A growing body of empirical evidence suggests that this ad hoc approach performs quite well. It has previously been argued that such a procedure works because it amounts to a sophisticated interpolation tool. We show that this is the case in a formal, asymptotic sense. In addition, we conduct some simulations which allow us to examine the importance of the sample size, the order of the polynomial, and the recalibration frequency in controlled settings. We also apply the ABS approach to daily S&P 100 index options to show that the procedure outperforms the Black-Scholes formula in valuing actual option prices out-of-sample.

Download the PDF at http://www.uh.edu/~jberkowi/ and the matlab files at http://www.bepress.com/snde/vol14/iss1/art4/.

Tags - black scholes , volatility

h = sqrt(kappa^2 + 2*sigma^2);

A = (2*h*exp((kappa+h)*T/2)/(2*h + (kappa+h)*(exp(T*h)-1)))^((2*kappa*alpha)/sigma^2);

B = 2*(exp(T*h)-1)/(2*h + (kappa+h)*(exp(T*h)-1));

P = A*exp(-B*r0); % bond price at 0

A = (2*h*exp((kappa+h)*T/2)/(2*h + (kappa+h)*(exp(T*h)-1)))^((2*kappa*alpha)/sigma^2);

B = 2*(exp(T*h)-1)/(2*h + (kappa+h)*(exp(T*h)-1));

P = A*exp(-B*r0); % bond price at 0

Tags - cir

Unfortunately, the authors don't release their codes for us to study, I tried to program according to that paper with theta scheme finite difference, where theta =0, 0.5, 1 refer to explicit, Crank-Nicolson, and implicit finite difference, respectively. Below is a runnable

% set parameter

N = 201; M = 200; s = 10; T = 1;

Tau = 0.1; %barrier window 20 days

sigma = 0.2; r = 0.05; K = 10;

Bar = 12; %barrier

bar = log(Bar);

% time-space grid

R = 3;

h = 2*R/(N+1);

k = T/M;

NoTau = floor(Tau/k);

x = linspace(-R,R,N+2)';

% compute finite difference matrix A

e = ones(N,1);

alphap = -sigma^2/2/h^2 +(sigma^2/2-r)/2/h;

alpham = -sigma^2/2/h^2 -(sigma^2/2-r)/2/h;

beta = sigma^2/h^2+r;

A = spdiags([alpham*e, beta*e, alphap*e], -1:1, N, N);

% compute matrices for the theta scheme

theta = 0.5;

B = speye(N,N) + theta*k*A;

C = speye(N,N) - (1-theta)*k*A;

% compute initial data

u = max(exp(x)-K,0);

u = repmat(u,1,NoTau+1);

inx = find(x>bar,1,'first');

u(:,1) = 0; %up and out when j=tau

f = zeros(N,1);

% start timestepping

for m = 1:M

lastu = u;

% compute right hand side

for j = 2:NoTau+1

f = C*[lastu(2:inx-2,NoTau+1);lastu(inx-1:end-1,j-1)]; %below using tau=0, above using tau=tau+1;

u(:,j) = zeros(N+2,1);

% solve the linear system

u(2:N+1,j) = B\f;

end

u(inx-1,2:NoTau) = u(inx-1,NoTau+1);%reset value at barrier point for parisian

lastu = u;

end

N = 201; M = 200; s = 10; T = 1;

Tau = 0.1; %barrier window 20 days

sigma = 0.2; r = 0.05; K = 10;

Bar = 12; %barrier

bar = log(Bar);

% time-space grid

R = 3;

h = 2*R/(N+1);

k = T/M;

NoTau = floor(Tau/k);

x = linspace(-R,R,N+2)';

% compute finite difference matrix A

e = ones(N,1);

alphap = -sigma^2/2/h^2 +(sigma^2/2-r)/2/h;

alpham = -sigma^2/2/h^2 -(sigma^2/2-r)/2/h;

beta = sigma^2/h^2+r;

A = spdiags([alpham*e, beta*e, alphap*e], -1:1, N, N);

% compute matrices for the theta scheme

theta = 0.5;

B = speye(N,N) + theta*k*A;

C = speye(N,N) - (1-theta)*k*A;

% compute initial data

u = max(exp(x)-K,0);

u = repmat(u,1,NoTau+1);

inx = find(x>bar,1,'first');

u(:,1) = 0; %up and out when j=tau

f = zeros(N,1);

% start timestepping

for m = 1:M

lastu = u;

% compute right hand side

for j = 2:NoTau+1

f = C*[lastu(2:inx-2,NoTau+1);lastu(inx-1:end-1,j-1)]; %below using tau=0, above using tau=tau+1;

u(:,j) = zeros(N+2,1);

% solve the linear system

u(2:N+1,j) = B\f;

end

u(inx-1,2:NoTau) = u(inx-1,NoTau+1);%reset value at barrier point for parisian

lastu = u;

end

Plots of the

Tags - parisian , finite-difference

Quotation

Matlab codes can be downloaded at http://www.rmi.nus.edu.sg/DuanJC/l, several other programming files can be also found at the page, for example, Co-integration option pricing model,

Tags - garch , volatility

My example in that post used simply Matlab embedded command "quadl", however, when the option becomes exotic, or the option has multiple observation times (for example, a Bermudan option), quadl can't be applied anymore as we have to match the nodes of these different observation times. (quadl does not allow users to specify how many points and how small interval an integral can have, correct me if i am wrong). Here is a good site on

Have fun at http://numericalmethods.eng.usf.edu/index.html. A minor defect is its gauss quadrature point is up to 10 only, no worry, should you are unhappy with that, another C/C++ code on Gauss quadrature can be found at http://www.holoborodko.com/pavel/?page_id=679.

Tags - quadrature , integral

For s short comparison, a simple QUAD code to price a vanilla European call option is as follows, please refer to the original paper for the meaning of symbols:

%using QUAD to calculate a vanila european call

x0 = log(s/k);

kappa = 2*(r-d)/vol^2-1;

dt = t;

ymax = 3;

A = 1/(sqrt(2*vol^2*pi*dt))*exp(-(kappa*x0/2)-(vol^2*kappa^2*dt/8)-r*dt);

q = quadl(@myquad,0,ymax,[],[],x0,dt,vol,kappa,k); %Matlab embedded quadrature

callprice = A*q;

function f = myquad(x,x0,dt,vol,kappa,k)

B = exp(-((x0-x).^2./(2*vol^2*dt))+kappa*x/2);

f = B.*max(exp(x)-1,0)*k;

x0 = log(s/k);

kappa = 2*(r-d)/vol^2-1;

dt = t;

ymax = 3;

A = 1/(sqrt(2*vol^2*pi*dt))*exp(-(kappa*x0/2)-(vol^2*kappa^2*dt/8)-r*dt);

q = quadl(@myquad,0,ymax,[],[],x0,dt,vol,kappa,k); %Matlab embedded quadrature

callprice = A*q;

function f = myquad(x,x0,dt,vol,kappa,k)

B = exp(-((x0-x).^2./(2*vol^2*dt))+kappa*x/2);

f = B.*max(exp(x)-1,0)*k;

Here I arbitrarily set ymax=3, which is enough for this simple example, the result for a European call option with strike price 9, stock price 10, volatility 20%, risk free rate 2%, dividend 1%, time to maturity 2 years is

The exotic option pricing is left for further experiment, have a nice evening.

PS: i was seriously drunken last weekend, my poor stomach.

Tags - quadrature , integral

Anyway, it is your turn to compare them. I write M files of several selected

pos=zeros(size(price,1),1);

[lead,lag]= movavg(price,x,y,'e');

lead = [zeros(x-1,1); lead]; %to avoid dimension mismatch

lag = [zeros(y-1,1); lag];

pos(lead>lag)=1;

pos(lag>lead)=-1;

[lead,lag]= movavg(price,x,y,'e');

lead = [zeros(x-1,1); lead]; %to avoid dimension mismatch

lag = [zeros(y-1,1); lag];

pos(lead>lag)=1;

pos(lag>lead)=-1;

m = size(price,1);

pos=zeros(m,1);

for i = 2:m

%put here the way to calculate variance C

UpperTrigger = price(i-1)+multiplier*sqrt(C);

LowerTrigger = price(i-1)-multiplier*sqrt(C);

if price(i)>=UpperTrigger pos(i) = 1;

elseif price(i)<=LowerTrigger pos(i) = -1;

end

end

pos=zeros(m,1);

for i = 2:m

%put here the way to calculate variance C

UpperTrigger = price(i-1)+multiplier*sqrt(C);

LowerTrigger = price(i-1)-multiplier*sqrt(C);

if price(i)>=UpperTrigger pos(i) = 1;

elseif price(i)<=LowerTrigger pos(i) = -1;

end

end

stosc = stochosc(highp, lowp, closep, kperiod, dperiod); %embedded Matlab function

m = size(highp,1);

pos=zeros(m,1);

inx1 = find(stosc(:,1)>=30);

inx2 = find(stosc(:,1)>=stosc(:,2));

pos(intersect(inx1,inx2)) = 1;

inx1 = find(stosc(:,1)<=80);

inx2 = find(stosc(:,1)<=stosc(:,2));

pos(intersect(inx1,inx2)) = -1;

m = size(highp,1);

pos=zeros(m,1);

inx1 = find(stosc(:,1)>=30);

inx2 = find(stosc(:,1)>=stosc(:,2));

pos(intersect(inx1,inx2)) = 1;

inx1 = find(stosc(:,1)<=80);

inx2 = find(stosc(:,1)<=stosc(:,2));

pos(intersect(inx1,inx2)) = -1;

%divergence index strategy, m is long momentum period, n is for short

longmom = tsmom(price,m);

shortmom = tsmom(price,n);

mm = size(price,1);

pos=zeros(mm,1);

DI = longmom.*shortmom./var(diff(price));

inx1 = find(DI<-8);

inx2 = find(longmom<0);

inx3 = find(longmom>0);

pos(intersect(inx1,inx2))=-1;

pos(intersect(inx1,inx3))=1;

longmom = tsmom(price,m);

shortmom = tsmom(price,n);

mm = size(price,1);

pos=zeros(mm,1);

DI = longmom.*shortmom./var(diff(price));

inx1 = find(DI<-8);

inx2 = find(longmom<0);

inx3 = find(longmom>0);

pos(intersect(inx1,inx2))=-1;

pos(intersect(inx1,inx3))=1;

% p - *price* data

% N - number of points to generate signal

pos=zeros(size(p,1),1);

macs = zeros(size(p,1),20);

for i=1:20

j=i*4;

[lead,lag]= movavg(p,i,j,'e');

lead = [nan(i-1,1); lead]; %to avoid dimension mismatch

lag = [nan(j-1,1); lag];

macs(lead>lag,i) = 5;

end

macssum = sum(macs,2);

macssum(1:80) = 50; %first 80 observations with zero position

pos(macssum>=N)=1;

pos(macssum<=(100-N))=-1;

% N - number of points to generate signal

pos=zeros(size(p,1),1);

macs = zeros(size(p,1),20);

for i=1:20

j=i*4;

[lead,lag]= movavg(p,i,j,'e');

lead = [nan(i-1,1); lead]; %to avoid dimension mismatch

lag = [nan(j-1,1); lag];

macs(lead>lag,i) = 5;

end

macssum = sum(macs,2);

macssum(1:80) = 50; %first 80 observations with zero position

pos(macssum>=N)=1;

pos(macssum<=(100-N))=-1;

There are other strategies left to you to backtest the effectiveness of technical analysis, for example, Kestner’s Moving Average System, Second Order Breakout, MACD Histogram Retracement, Normalized Envelope Indicator, etc. Have fun.

Tags - trading , strategy

Click here to watch it, including:

Quotation

• building a robust and generic back-testing framework

• building and testing nonlinear trading models

• using parallel computing to improve efficiency

• deploying models to work with third party software

• building and testing nonlinear trading models

• using parallel computing to improve efficiency

• deploying models to work with third party software

Tags - matlab , trading , algo

1, Writing Fast MATLAB Code, PDF doc;

2, Maximize Matlab performance, http://www.ece.northwestern.edu/local-apps/matlabhelp/techdoc/matlab_prog/ch7_per6.html;

3, MATLAB array manipulation tips and tricks, http://home.online.no/~pjacklam/matlab/doc/mtt/index.html;

4, matlab tips and tricks, http://www.ee.columbia.edu/~marios/matlab/matlab_tricks.html;

5, Accelerating Matlab, http://research.microsoft.com/en-us/um/people/minka/software/matlab.html

Many other Matlab accelerating skills can be found on the webpage 4. Be ready to take off.

Tags - matlab

2: Neural network;

3: Multiple regression model;

4: A very simple model based on a simple moving average.

What interests me is

(kernel or Gram) matrix. For our aims it is natural to think of the price of a stock as being some function over time. Generally, a Gaussian processes can be cosidered to be defining a distrubution over functions with the inference step occurring directly in the space of functions. Thus, by using

Should you are interested, here is a book Gaussian Processes for Machine Learning to be freely downloaded, accompanying Matlab package is also available at the website.

http://www.gaussianprocess.org/gpml/

Tags - regression

Quantitative Trading: How to Build Your Own Algorithmic Trading Business is a great book I bought several months ago but had few chances to read until recently. It is an easy-to-understand book on implementing quantitative (or algorithmic) trading for independent traders. This is the area fascinating me from right the beginning of learning quantitative finance. Therefore I have started to play

Recorded Webinar: Algorithmic Trading with MATLAB for Financial Applications: http://www.mathworks.com/company/events/webinars/wbnr30376.html?id=30376&p1=50647&p2=50649

Interactive Brokers via Matlab: http://www.maxdama.com/2008/12/interactive-brokers-via-matlab.html

There you could find and download useful resources and API code on

Tags - matlab , trading

Quotation

A new methodology is proposed to estimate theoretical prices of financial contingent claims whose values are dependent on some other underlying financial assets. In the literature, the preferred choice of estimator is usually maximum likelihood (ML). ML has strong asymptotic justification but is not necessarily the best method in finite samples.

This paper proposes a simulation-based method. When it is used in connection with ML, it can improve the finite-sample performance of the ML estimator while maintaining its good asymptotic properties. The method is implemented and evaluated here in the Black-Scholes option pricing model and in the Vasicek bond and bond option pricing

model. It is especially favored when the bias in ML is large due to strong persistence in the data or strong nonlinearity in pricing functions. Monte Carlo studies show that the proposed procedures achieve bias reductions over ML estimation in pricing contingent claims when ML is biased. The bias reductions are sometimes accompanied by reductions in variance. Empirical applications to U.S. Treasury bills highlight the differences between the bond prices implied by the simulation-based approach and those delivered by ML. Some consequences for the statistical testing of contingent-claim pricing models are discussed.

This paper proposes a simulation-based method. When it is used in connection with ML, it can improve the finite-sample performance of the ML estimator while maintaining its good asymptotic properties. The method is implemented and evaluated here in the Black-Scholes option pricing model and in the Vasicek bond and bond option pricing

model. It is especially favored when the bias in ML is large due to strong persistence in the data or strong nonlinearity in pricing functions. Monte Carlo studies show that the proposed procedures achieve bias reductions over ML estimation in pricing contingent claims when ML is biased. The bias reductions are sometimes accompanied by reductions in variance. Empirical applications to U.S. Treasury bills highlight the differences between the bond prices implied by the simulation-based approach and those delivered by ML. Some consequences for the statistical testing of contingent-claim pricing models are discussed.

Download the paper and accompanying matlab code at http://www.mysmu.edu/faculty/yujun/research.html

Tags - simulation , mle

for forecasting. Wikipedia has a detailed explanation on it.

Unfortunately I have not used it except once I tried the built-in VAR function in Eviews over 5 years ago, when one of my classmates did a seminar presentation on it. Sorry I couldn't find useful VBA code, what i do get is a sample spreadsheet showing the VAR Series Analysis & Results but it seems the author intentionly hides the macro code, http://www.afpc.tamu.edu/courses/622/files/lecturedemos/Lecture%2007%20Vector%20Autoregression.xls.

If you are happy with Matlab, here is a

Tags - vector-autoregression , var

Audio - Astronomy - BiomedicalInformatics - Chemometrics - Chaos - Chemistry - Coding - Control - Communications - Engineering - Data Mining - Excel - FEM - Fuzzy - Finance - GAs - Graph - Graphics - Images - ICA - Kernel - Markov - Medical - MIDI - Misc. - MPI - NNets - Oceanography - Optimization - Plot - Signal Processing - Optimization - Statistics - SVM - Web - etc ...

Recommended

Kernel Density Estimation Toolbox

http://ssg.mit.edu/~ihler/code/

BOOTSTRAP MATLAB TOOLBOX

http://www.csp.curtin.edu.au/downloads/bootstrap_toolbox.html

CompEcon Toolbox for Matlab

http://www4.ncsu.edu/~pfackler/compecon/toolbox.html

Random Neural Networks

http://www.cs.ucf.edu/~ahossam/rnnsimv2/

Logistic regression

http://www.spatial-econometrics.com/

ARfit: A Matlab package for the estimation of parameters and eigenmodes of multivariate autoregressive models

http://www.gps.caltech.edu/~tapio/arfit/

Time Series Analysis

http://www.dpmi.tu-graz.ac.at/~schloegl/matlab/tsa/

Interested ppl may download more at http://www.tech.plym.ac.uk/spmc/links/matlab/matlab_toolbox.html

Tags - matlab , toolbox

steep.m : Steepest Descent

gaussn.m : Damped Gauss-Newton

bfgswopt.m : BFGS, low storage

Polynomial line search routines: polyline.m , polymod.m

Numerical Derivatives: diffhess.m : Difference Hessian,

requires dirdero.m : directional derivative, as do several other codes

ntrust.m : Newton's Method with Simple Dogleg

levmar.m : Levenberg-Marquardt for nonlinear least squares

cgtrust.m : Steihaug CG-dogleg

gradproj.m : Gradient Projection Method

projbfgs.m: Projected BFGS code

imfil.m : Implicit Filtering

nelder.m : Nelder-Mead

simpgrad.m : Simplex Gradient, used in implicit filtering and Nelder-Mead codes

hooke.m : Hooke-Jeeves code

mds.m : Multidirectional Search code

Check http://www.siam.org/books/kelley/fr18/matlabcode.php for detail.

Tags - optimization

An interesting paper "A Data-Driven Optimization Heuristic for Downside Risk Minimization" demonstrates how to apply

Tags - heuristic , optimization

Here is a paper on the implementation of

Tags - binomial , american

But how many ways are you able to

here are the paper and matlab codes, you might feel in the end

Enjoy.

Tags - binomial

Download the

Download TradingSolutions

Tags - neural-network

Quotation

We propose a unified, fully general methodology to analyze and act on **diversification** in any environment, including long-short trades in highly correlated markets. First, we build the **diversification** distribution, i.e. the distribution of the uncorrelated bets in the portfolio that are consistent with the portfolio constraints. Next, we summarize the wealth of information provided by the **diversification** distribution into one single diversification index, the effective number of bets, based on the entropy of the **diversification distribution**. Then, we introduce the mean-diversification efficient frontier, a diversification approach to portfolio optimization. Finally, we describe how to perform mean-diversification optimization in practice in the presence of transaction and market impact costs, by only trading a few optimally chosen securities.

Click for codes http://www.mathworks.com/matlabcentral/fileexchange/23271

Tags - diversification