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];

1x4 struct array with fields:
    symbol
    desc
    lastTrade
    lastTradeTime
    dividendDate
    open
    lastClose

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

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

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

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] ======  

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.

%% 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)

% 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

a = rnorm(10,1);

b = rnorm(5,5);
c = rnorm(3,3)...
    + rnorm(3,3);

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

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

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];

% 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));

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);

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

% 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

%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;

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;

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

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;

%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;

% 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;