Codes collection for my formal projects in university

最新推荐文章于 2024-06-14 00:39:34 发布

weixin_30594001

最新推荐文章于 2024-06-14 00:39:34 发布

阅读量844

点赞数

文章标签： matlab 人工智能 c/c++

原文链接：http://www.cnblogs.com/jessewei/p/5035769.html

版权

In fact, these codes mean much more than a collection as they witnessed my growth.

There is a better version, a delectable PDF brochure, which you can download in

http://yun.baidu.com/share/link?shareid=4213336497&uk=539384257

I. Papers

Thanks to professors and teammates, I got the chances to train and exhibit my ability in modeling and programming. Here are two papers' codes.

A. Is China Stock Market Efficient

1. Introduction:

China stock market had last for 20 years. But we still wondered whether the market was rational. So Siyu Ye and I cooperated to analyze the data and published the paper, Is China Stock Market Efficient?—Basing on Improved BSJ Model. We examined whether China stock prices complied with the principles of CAPM and found the answer was NO.

2. Codes:

%20个组合，每个组合5支股票

%以regress的方式算出beta值,

%第一期涉及的数据为1----t1+1;第二期涉及到得数据为t1+1----t1+1+t2+2;第三期涉及到的数据为 t1+t2+2----end

% x为股票价格矩阵；rf为无风险利率列向量；index为市场指数，如上证综指

%构造股票超额收益率

X=diff(log(x))-repmat(rf(1:length(x)-1),1,size(x,2)); %差分减少了一个数据单位

%构造市场超额收益率

Index=diff(log(index))-rf(1:end-1);

H=[];%对时间序列进行平稳性检验，0为接受原假设，时间序列为非平稳序列；1为拒绝原假设

for i=1:100

h= pptest(X(:,i),'model','ar');

H=[H,h];

end

HI= pptest(Index,'model','ar');

%第一个时期开始

t1=ceil(length(X)/3);

B=[];%参数值

BINT0=[];%参数0（即常数项）的95%置信区间

BINT1=[];%参数1（即系数）的95%置信区间

Sr2=[];%R-square值

Sf=[];%f检验量

Sfp=[];%f检验量对应的P值

Std=[];%标准差

Resid=[];%残差

T=[];

TP=[];

JB=[];%jb检验结果 0为接受原假设，即95%把握为正态分布；1为不接受原假设

White=[];%white检验检验结果 0为接受原假设，即残差为同方差；1为不接受，即为异方差

TR2=[];%记录white检验中的TR2的值

DW=[];%自相关检验结果 0为接受原假设，即残差没有自相关；1为不接受

DWP=[];

%构造时期1的市场超额收益率和股票超额收益率

Index1=Index(1:t1,1);

X1=X(1:t1,:);

a=ones(length(Index1),1);%构造值为1的列向量

Index1=[a Index1];%构造自变量矩阵

for i=1:size(x,2) %股票有100支

xx=X1(:,i);

[b,bint,r,rint,statsl]=regress(xx,Index1);

%将股票超额收益率对市场超额收益进行回归

B=[B b];

BINT0=[BINT0 bint(:,1)];

BINT1=[BINT1 bint(:,2)];

Resid=[Resid r];

Sr2=[Sr2 statsl(:,1)];

Sf=[Sf statsl(:,2)];

Sfp=[Sfp statsl(:,3)];

Std=[Std statsl(:,4)];

h=jbtest(r); %JB检验

JB=[JB h];

qq=regstats(xx,Index1(:,2),'linear','tstat');

T=[T,qq.tstat.t(2)];

TP=[TP,qq.tstat.pval(2)];

qq=regstats(r,Index1(:,2),'linear','dwstat');%自相关检验

dwp=qq.dwstat.pval;

DWP=[DWP dwp];

dwp(dwp>0.05)=2;

dwp(dwp<0.05)=1;

dwp(dwp==2)=0;

DW=[DW dwp];

u=r.^2; %异方差white检验

[b,bint,r,rint,statsl]=regress(u,[Index1,Index1(:,2).^2]);

t=length(u);

tr2=t*statsl(1);

TR2=[TR2 tr2];

tr2(tr2<=chi2inv(0.95,2))=0;

tr2(tr2>chi2inv(0.95,2))=1;

White=[White tr2];

end

beta1=B(2,:); %选择第二行，即非常数项

%时期1结束

t2=t1;

[c,pos]=sort(beta1); %c为升序的beta值，pos为升序后的排序标记

%分成10个组合

s=zeros(length(x),20); %构造10个值为0的组合

for j=1:20 %加总以构造10个组合

for i=1:5 %10为每组的股票个数

s(:,j)=s(:,j)+x(:,pos(i+5*(j-1)));

end

s2=s(t1+1:t1+2+t2,:); %组合在第二期的价格变化

X2=diff(log(s2))-repmat(rf(t1+1:t1+t2+1),1,size(s2,2)); %差分减少了一个数据单位 X2是股票组合的超额收益率

Index2=Index(t1+1:t1+1+t2);

a=ones(length(Index2),1);%构造值为1的列向量

Index2=[a Index2];%构造自变量矩阵

pB=[];%参数值

pBINT0=[];%参数0（即常数项）的95%置信区间

pBINT1=[];%参数1（即系数）的95%置信区间

pSr2=[];%R-square值

pSf=[];%f检验量

pSfp=[];%f检验量对应的P值

pStd=[];%标准差

pResid=[];%残差

pT=[];

pTP=[];

pJB=[];

pWhite=[];

pTR2=[];

pDW=[];

pDWP=[];

for i=1:size(s2,2) %股票有100支

xx=X2(:,i);

[b,bint,r,rint,statsl]=regress(xx,Index2);

%将股票超额收益率对市场超额收益进行回归

pB=[pB b];

pBINT0=[pBINT0 bint(:,1)];

pBINT1=[pBINT1 bint(:,2)];

pResid=[pResid r];

pSr2=[pSr2 statsl(:,1)];

pSf=[pSf statsl(:,2)];

pSfp=[pSfp statsl(:,3)];

pStd=[pStd statsl(:,4)];

h=jbtest(r);

pJB=[pJB h];

qq=regstats(xx,Index2(:,2),'linear','tstat');

pT=[pT,qq.tstat.t(2)];

pTP=[pTP,qq.tstat.pval(2)];

qq=regstats(r,Index2(:,2),'linear','dwstat');

pdwp=qq.dwstat.pval;

pDWP=[pDWP, pdwp];

pdwp(pdwp>0.05)=2;

pdwp(pdwp<0.05)=1;

pdwp(pdwp==2)=0;

pDW=[pDW pdwp];

u=r.^2; %异方差white检验

[b,bint,r,rint,statsl]=regress(u,[Index2,Index2(:,2).^2]);

t=length(u);

ptr2=t*statsl(1);

pTR2=[pTR2 ptr2];

ptr2(ptr2<=chi2inv(0.95,2))=0;

ptr2(ptr2>chi2inv(0.95,2))=1;

pWhite=[pWhite ptr2];

end

beta2=pB(2,:); %beta2为10个组合的beta值

%到此步骤二结束

%第三期开始

s3=s(t1+2+t2:end,:); %组合在第3期的价格变化

portfolio=diff(log(s3)); %portfolio是股票组合的收益率

P=mean(portfolio); %组合的平均收益率，其对应的beta值为beta2

BB=[ones(length(beta2),1),beta2']; %构造自变量矩阵

[b,bint,r,rint,statsl]=regress(P',BB); %重要信息包含在statsl里

jb=jbtest(r);

qq=regstats(P',BB(:,2),'linear','tstat'); %t值

t=qq.tstat.t(2);

tp=qq.tstat.pval(2);

qq=regstats(r,BB(:,2),'linear','dwstat'); %自相关检验

dwp=qq.dwstat.pval;

dwp(dwp>0.05)=2;

dwp(dwp<0.05)=1;

dwp(dwp==2)=0;

u=r.^2; %异方差white检验

[fB,fBint,fR,fRint,fStatsl]=regress(u,[BB,BB(:,2).^2]);

tr2=length(u)*fStatsl(1);

tr2(tr2<=chi2inv(0.95,2))=0; %最后 tr2包含着最后的回归的异方差信息

tr2(tr2>chi2inv(0.95,2))=1;

B. Make money by a dynamic method

1. Introduction:

Last semester, Zhang Yumin and I cooperated to write a paper, the Best Portfolio under the Dynamic Risk Measurement of CVAR Model—Basinng on the Psychological Frame, in which I combined the dynamic optimum algorithm and the CVAR risk measurement in order to find some specific way to minimize the risk.

2. Codes:

a) Script 1:

function v =SubObjF2(k,x,u) %阶段k的指标函数

%q=wA历史模拟法

%S是244行4列的矩阵

load('s.mat');

retmin=[0.00220811 0.002253818 0.002619162 0.00358249

];

port=S([1:61]+(k-1)*61,:); %61也可以改<-----这个地方如果资产改的话这个代码要改

port=(x+u)*S';%得出一整行的数据其中 x+u 是关键

portret=sum(port);

if portret>retmin(k)

a=sort(port);

t=length(a);

num=ceil(0.05*t);%求出前5%小的数在a中的位置

portvar=a(num);

portcvar=mean(min(port,portvar)-portvar)*1/(1-0.95)+portvar;%性能指标

v=-portcvar; %函数求最小

else

v=inf;

end

b) Script 2:

function [p_opt, fval] = dynprog2(x,DecisF2, SubObjF2, TransF2, ObjFun)

% 动态规划逆序算法

% x是状态变量，一列代表一个阶段状态；

% M-函数DecisFun(k,x)由阶段k的状态变量x求出相应的允许决策变量；

% M-函数SubObjFun(k,x,u)是阶段指标函数;

% M-函数TransFun(k,x,u)是状态转移函数，其中x是阶段k的某状态变量，u是相应的决策变量；

% M-函数ObjFun(v,f)是第k阶段至最后阶段指标函数，当ObjFun(v,f) = v+f时，输入的ObjFun可以省略；

% fval是一个列向量，各元素分别表示p_opt各最优策略组对应始端状态x的最优函数值

k =size(x,3); %k个时间阶段

%x_isnan = ~isnan(x); %非空状态坐标 %本题应用中并不存在空状态

t_vubm = inf*ones(size(x,1),k); %性能指标中间矩阵

f_opt = nan*ones(size(x,1),k); %总性能指标矩阵

d_opt =nan*ones(size(x)); %每步决策矩阵并存储最优决策

%计算终端（最后一阶段）相关值

%tmp1 = find(x_isnan(:,k)); %最后一步状态向量

tmp2 = size(x,1); %状态数量

for i = 1:tmp2

t_vub = inf; %%用于比较各性能指标，并存储最优性能指标值

u = feval(DecisF2, k,x(i,:,k)); %决策变量

tmp3 = size(u,1); %策略数量

for j = 1:tmp3 %求出当前状态下所有决策的最小性能指标

tmp = feval(SubObjF2,k,x(i,:,k),u(j,:));

if tmp <= t_vub

f_opt(i,k) = tmp;

d_opt(i,:,k) = u(j,:);

t_vub =tmp;

end

for ii = k-1:-1:2 %%%%%%%原本2是1

%tmp10 = find(x_isnan(:,ii));

tmp20 = size(x,1);

for i = 1:tmp20 %求出当前状态下所有可能的决策

u = feval(DecisF2,ii,x(i,:,ii));

tmp30 = size(u,1);

for j = 1:tmp30 %求出当前状态下所有可能决策的最小性能指标

tmp00 = feval(SubObjF2, ii, x(i,:,ii),u(j,:)); %当步性能指标

tmp40 = feval(TransF2,ii,x(i,:,ii),u(j,:)); %下一状态向量

%弃用unique

%xtmp=[tmp40;x(:,:,1)];

%[C,ia,ic]=unique(xtmp,'rows');

%ic(1)=[];

%position=find(ic==1); %找出下一状态在x矩阵的位置

tmp4=repmat(tmp40,size(x,1),1);

tmp5=x(:,:,1)-tmp4;

tmp5(abs(tmp5)<0.0001)=0;

b=logical(tmp5);

sumb=sum(b,2);

position=find(sumb==0); %找出下一状态在x矩阵的位置

tmp50 = x(position,:,ii+1);

if ~isempty(tmp50)

if nargin < 5

tmp00 =tmp00 + f_opt(position,ii+1);

else

tmp00 = feval(ObjFun,tmp00,f_opt(position,ii+1));

end %%当前状态的性能指标

if tmp00 <= t_vubm(i,ii)

f_opt(i,ii) = tmp00;

d_opt(i,:,ii) = u(j,:);

t_vubm(i,ii) = tmp00;

end

ii = 1 ;

xx=[0.5,0.2,0.1,0.2];%%%%%%%%%%%%%%%%%%%初始状态

%tmp10 = find(x_isnan(:,ii));

%求出当前状态下所有可能的决策

u = feval(DecisF2,ii,xx);

tmp30 = size(u,1);

for j = 1:tmp30 %求出当前状态下所有可能决策的最小性能指标

tmp00 = feval(SubObjF2, ii, x(1,:,ii),u(j,:)); %当步性能指标

tmp40 = feval(TransF2,ii,x(1,:,ii),u(j,:)); %下一状态向量

%弃用unique

%xtmp=[tmp40;x(:,:,1)];

%[C,ia,ic]=unique(xtmp,'rows');

%ic(1)=[];

%position=find(ic==1); %找出下一状态在x矩阵的位置

tmp4=repmat(tmp40,size(x,1),1);

tmp5=x(:,:,1)-tmp4;

tmp5(abs(tmp5)<0.0001)=0;

b=logical(tmp5);

sumb=sum(b,2);

position=find(sumb==0); %找出下一状态在x矩阵的位置

tmp50 = x(position,:,ii+1);

if ~isempty(tmp50)

if nargin < 5

tmp00 =tmp00 + f_opt(position,ii+1);

else

tmp00 = feval(ObjFun,tmp00,f_opt(position,ii+1));

end %%当前状态的性能指标

if tmp00 <= t_vubm(1,ii)

f_opt(1,ii) = tmp00;

d_opt(1,:,ii) = u(j,:);

t_vubm(1,ii) = tmp00;

end

fval = f_opt;

p_opt=d_opt;

%fval = fval(find(~isnan(fval)),1);

%记录最优决策、最优轨迹线和相应指标函数值。

c) Script 3:

function u = DecisF2(k,x) %在阶段k由状态变量x的值求出其相应的决策变量所有的取值

temp=size(x);

uup=ones(temp);

udown=uup;

bound=0.2;

step=0.1;

high=1-x; %1是单资产的持有上限，high是上升幅度上限

uup(high<bound)=high(high<bound);

uup(high>=bound)=bound; %策略增量

low=x-0; %0是单资产的持有下限 low是下降幅度上限

udown(low<bound)=low(low<bound);

udown(low>=bound)=bound; %策略减少量

[x1 x2 x3 x4]=ndgrid(-udown(1):step:uup(1),-udown(2):step:uup(2),-udown(3):step:uup(3),-udown(4):step:uup(4));

thesum=x1+x2+x3+x4;

ix=find(abs(thesum)<eps);

u=[x1(ix),x2(ix),x3(ix),x4(ix)];

d) Script 4:

function y = TransF2(k,x,u) %状态转移方程

y = x+u;

e) Script 5:

clear

%求解所有状态空间

tic

n=100; %单个资产上限

gap=10; %状态间隔

[x1 ,x2 ,x3 ,x4]=ndgrid(0:gap:n,0:gap:n,0:gap:n,0:gap:n);

thesum=x1+x2+x3+x4;

ix=find(abs(thesum-100)<eps);

answ=[x1(ix),x2(ix),x3(ix),x4(ix)];

%总的阶段数为T

T=4; %

x=repmat(answ,[1,1,T]);%x是状态总空间每一列为一种权重分配，即一个向量%%%%x=answ;%

x=x/100;

%原问题中的x是一个9行3列的矩阵，每一列代表一个阶段的状态空间，而我们这里就省略了，直接引用状态总空间，故而在主函数部分需要较大的改动

[p,f] = dynprog2(x,'DecisF2','SubObjF2','TransF2')

a=[0.5,0.2,0.1,0.2];%%%%%%%%%初始状态

u=p(1,:,1);

Posi=[];

for k=2:T

a(k,:)=a(k-1,:)+u(k-1,:);

tmp4=repmat(a(k,:),size(x,1),1);

tmp5=x(:,:,1)-tmp4;

tmp5(abs(tmp5)<0.0001)=0;

b=logical(tmp5);

sumb=sum(b,2);

position=find(sumb==0); %阶段k时最佳状态所在位置

Posi=[Posi;position];

u(k,:)=p(position,:,k);

end

a(5,:)=a(4,:)+u(4,:);

toc

%得出a是总状态，u是总策略

II. Mathematical Modeling Contest

No wonder I got extensive progress during the preparation period for the Mathematic Modeling Contest, for not only the books I read but the experience I obtained from modeling.

I attended totally 2 such contests and here are the codes.

A. Fight against Ebola

1. Introduction:

Ebola virus had spread many countries and threatened people there. In February 2015, we participated in the North America Mathematic Modeling Contest and Ebola was one of the questions. So we designed a model to simulate the spread process of virus by using Monte Carlo Simulation. Basing on the type of the spread, we continued to apply K-Mean method to calculate a fastest distribution that delivered the medicine to the dire places.

2. Codes:

a) Script 1:

%核心模型的修缮

%已经加入随机项

%进行蒙特卡洛模拟法

%一定数量的药物所带来的，模拟非最优方式下的结果，设置固定的随机方式

%不以比例，而以数量进行计算

clear;close all;

%这个拟合不以比例进行

%三个国家按顺序分别是几内亚利比里亚塞拉利昂

P=[];

N=[11.75;4.294;6.092]*10^6;

load('shumo.mat');%载入mat

load('Time.mat');%载入Time

% Time=(1:size(mat,1));

for i=1:3

perc(:,i)=mat(:,i);%患者数量

end

y=log(perc);

index=y~=-inf;

for i=1:3

xtemp=Time(index(:,i));

ytemp=y(:,i);

ytemp=ytemp(index(:,i));

p=polyfit(xtemp,ytemp,1);

P=[P;p];%拟合出beta*N值，为P的第一列

end

%用以连接下一个m文件

beta=P(:,1)./N;%三个国家的beta值，疾病扩散系数

I0=perc(end,:)';%引进药物以后的初始数

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

PERCT=[];

PSN=[];TT=[];

%我们需要的变量，在循环前初始化

for l=1:100

n=3;%假设的三个不同的城市，视为质点

t=30;%时间长度，天

positions=rand(n,2);%城市的位置

%plot(positions(:,1),positions(:,2),'o');

I=I0;

%或者进行修改 I0N=[2975-1944;8745-3746;10740-2376];I=I0N./N;

S=N-I;

gamma=0.0000001;

%p=zeros(n,1);%不同城市的药物分配量

R=beta./gamma.*S;%交换数，一个病人的康复期内与导致的染病人数的比值

%该值小于1，说明疾病能够得到遏制

%perct=rand(n,1);%进行随机

%perct=N.*N;%停止随机

perct=N.^3;

PERCT=[PERCT,perct];

prio=(R.*I).*perct;%作为一个尝试的分配标准 %perct=N.^3 2.210735

pmax=1500;%设置一个药厂的最高产量。同时药厂还需要一个初始最低产量以保证疾病的不扩散

%随着药量的改变，pSN的变化回事如何，T的变化有事如何

p=pmax*prio/nansum(prio);%不同城市的药物分配量

P=[];

P=[P,p];

GAMMA=[];

%添加随机项e，以二项分布模拟，取值与每座城市染病人数有关

mu=0.1; %产生随机项的重要系数

E=[];

for i=1:t

for j=1:n %第j个城市的随机项

r=binornd(floor(mu*I(j,i)),0.5);

r=ceil(r);%e必然大于等于0

temp=sign(randn);

e(j)=r*temp;%e随机正负

end

E=[E,e'];

e=0;%关闭随机项

S(:,i+1)=S(:,i)-beta.*S(:,i).*I(:,i)+e';

S(:,i+1)=-min(0,-S(:,i+1));%S不低于0

I(:,i+1)=I(:,i)+beta.*S(:,i).*I(:,i)-p-e';

I(:,i+1)=-min(0,-I(:,i+1));%I不低于0

gamma=p./(N.*I(:,i+1));

R=beta./gamma.*S(:,i+1);

prio=R.*I(:,i+1).*perct;

p=pmax*prio/nansum(prio);%不同城市的药物分配量

GAMMA=[GAMMA,gamma];

P=[P,p];

end

P(isnan(P))=0;

a=(I==0);

s=sum(a);

T=find(s==n, 1 );%第T天我们所有城市地患病人消失，即瘟疫结束，费时T-1天

if isempty(T)

display('瘟疫在给定时间内不停止');

T=inf;

else

%接下来，通过评估第T天的数据，我们就能确定我们的方案是否是最优的

pSN=sum(S(:,T));%瘟疫结束后，总的未感染的人

pIN=N-pSN;

pP=P(:,1:T-1);

end

PSN=[PSN;pSN];

TT=[TT;T];

end

[TT,PSN];

mean(PSN)

b) Script 2:

%数据准备与初始化

%初级聚类点可以作为次级聚类点

clear;close;

x=randn(30,2);

n=length(x);%感染地数量

medicine=abs(1000+400*randn(n,1));%每个城市在整个过程中总的药物的分配量

medic=diag(medicine);

y=pdist(x);

s=squareform(y);%各点之间的欧式距离矩阵

sini=s;%记录下最初的数据

vsub=1;%定义支路运输速度

vprim=1;%定义干路运输速度

%寻找初级聚类中心

DDtotal=[];

for prim=1:n

medic=1;%将药量考虑进模型的开关

s=sini*medic;

%s(:,prim)=[];

%s(prim,:)=[];

%寻找次级聚类中心

Dtotal=[];

for pa=1:n-1 %a在b前，同时又要减去一个初级聚类中心，故循环到n-2

for pb=pa+1:n %我们选取(去掉第prim个城市后的地图中)第pa个城市（城市a）和第pb个城市（城市b）为次级聚类中心

dprima=sini(prim,pa);

dprimb=sini(prim,pb);

dprim=dprima+dprimb;%干路总距离

a=s(pa,:);%第p1个城市到其他城市的距离

b=s(pb,:);%第p2个城市的距离

comp=(a<=b);%logical变量

da=a(comp);

db=b(~comp);

dsuma=sum(da);

dsumb=sum(db);

dmeana=dsuma/max((sum(comp)-1),1);%使用max，使得作为除数的待分配城市数不至于为0

dmeanb=dsumb/max((sum(~comp)-1),1);

dtotal=max(max(da/vsub+dprima/vprim),max(db/vsub+dprimb/vprim));%从时间最优的角度出来的标准

%另外，可以添加一个从道路运输负担出来的标准

Dtotal=[Dtotal;dtotal];%记录下在选定prim-delivery点时，不同sub-delivery点所得到的标准

%通过比较该列向量，可以得出某一最佳的sub-delivery的点

end

%这里不同的列记录下不同的prim点会带来的最优点的集合

DDtotal=[DDtotal,Dtotal];

end

m=min(min(DDtotal));

[row,col]=find(DDtotal==m);%第col个城市作为prim-delivery是最好的

col=col(1);

fd=DDtotal(:,col);

fdsquare=squareform(fd);

[frow,fcol]=find(fdsquare==m,1,'last');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%---------------------------

%找到最终聚类comp

%可以的出各种关于最佳聚类的数据！！！

prim=col;

medic=1;%将药量考虑进模型的开关

s=sini*medic;

%寻找次级聚类中心

Dtotal=[];

pa=frow; %a在b前，同时又要减去一个初级聚类中心，故循环到n-2

pb=fcol; %我们选取(去掉第prim个城市后的地图中)第pa个城市（城市a）和第pb个城市（城市b）为次级聚类中心

dprima=sini(prim,pa);

dprimb=sini(prim,pb);

dprim=dprima+dprimb;%干路总距离

a=s(pa,:);%第p1个城市到其他城市的距离

b=s(pb,:);%第p2个城市的距离

comp=(a<=b);%logical变量去掉初级聚类城市后的聚类参数

da=a(comp);

db=b(~comp);

dsuma=sum(da);

dsumb=sum(db);

dmeana=dsuma/max((sum(comp)-1),1);%使用max，使得作为除数的待分配城市数不至于为0

dmeanb=dsumb/max((sum(~comp)-1),1);

dtotal=max(max(da/vsub+dprima/vprim),max(db/vsub+dprimb/vprim));%从时间最优的角度出来的标准

lll=1:n;

citya=lll(comp);%以第numofa个城市为中心的a类城市群

cityb=lll(~comp);%以第numofb个城市为中心的b类城市群

cityprim=col;%prime-delivery

%%%%%%%%%%%%%%%%%%%%%%%----------------------------

xx=x(:,1);

yy=x(:,2);

plot(xx,yy,'k.'); %画出每个城市在地图上的点

hold on;

for i=1:n %标号

c=num2str(i);

c=[' ',c];

text(xx(i),yy(i),c,'fontsize',18);

end

plot(xx(citya),yy(citya),'bo','markersize',8,'markerfacecolor','b');

plot(xx(cityb),yy(cityb),'ms','markersize',8,'markerfacecolor','m');

plot(xx(col),yy(col),'r.','markersize',10);

%初级聚类中心

plot(xx(col),yy(col),'kh','markersize',25);

%次级聚类中心

plot(xx(frow),yy(frow),'+','markersize',30,'markeredgecolor','b');

plot(xx(fcol),yy(fcol),'+','markersize',30,'markeredgecolor','m');

B. Find your location even without any modern devices.

1. Introduction:

This is the algorithm for the Chinese University Students' Mathematic Modeling Contest. Our mission was to calculate the location of a vertical stick with its shadow movemonts being recorted every 20 minutes. I solve it with the PSO-AG method, which is originally combine the particle swarm optimization with genetic algorithm. And in the further exploration, we even provided the date for those these data. So you could still know where and when it is even you were lost in a desolated island. Cool.

(I leave out some codes that has the similar function with following codes. If you want to overview the whole codes, please open the attachments)

2. Codes:

a) Script 1:

% lati 纬度

% decli 赤纬

% longi 经度

% hour_angle 时角

% h 高度角

% A 方位角

% sin(h)=sin(alti)*sin(decli)+cos(alti)*cos(decli)*cos(hour_angle)

% cos(A)=( sin(decli)*sin(alti)-cos(decli)*cos(alti)*cos(hour_angle)

% )/sin(h)

lati=(39+54/60+26/60/60);

longi=(116+23/60+29/60/60); %这里暂时是角度制

lati=lati/180*pi;

%n为天数从1月1号开始计算

N=datenum(2015,10,22)-datenum(2015,1,1)+1;

% decli=( 23.45*sin(360*(284+N)/365 /180*pi) )/180*pi;

b=2*pi*N/365;

decli=( 0.006918-0.399912*cos(b)+0.070257*sin(b)-0.006758*cos(2*b)+0.000907*sin(2*b)-0.002697*cos(3*b)+0.00148*sin(3*b) );

%t为北京时间 tt为当地时间即真太阳时

t=9: 1/60 :15;

tt=t+(longi-120)*4/60;

hour_angle=(tt-12)*15/180*pi;

sin_h=sin(lati)*sin(decli)+cos(lati)*cos(decli)*cos(hour_angle);

cos_A=( sin_h * sin(lati) - sin(decli) )./ (sqrt(1-sin_h.^2)*cos(lati));

% cos_A=( sin(decli)*sin(lati)-cos(decli)*cos(lati)*cos(hour_angle) )./sin_h;

L=3;

r=L*( sqrt(1-sin_h.^2)./sin_h ); %部分点乘被省略，注意添加

% x=r.*cos_A;

% y=r.*sqrt(1-cos_A.^2);

%其中delta的单位为度(deg)；pi=3.1415926为圆周率；b(deg)=2*pi*(N-1)/365，单位为度(deg)；

% 0.39795*cos(0.98563*(N-173))

plot(t,r);

legend('影子长度随时间变化示意图')

xlabel('时间 /时')

ylabel('影子长度 /米')

b) Script 2:

E0=0.000000001; %允许误差

MaxNum=200; %粒子最大迭代次数

narvs=3; %%%%%%% %目标函数的自变量个数

particlesize=50; %粒子群规模

c1=2; %每个粒子的个体学习因子，也称为加速常数

c2=2; %每个粒子的社会学习因子，也称为加速常数

w=0.75; %惯性因子

vmax=0.8; %粒子的最大飞翔速度

nn=30; % 用于遗传算法的初始种群数量

AGx=ones(nn,narvs);

AGfitness=ones(nn,1);

for q=1:nn

% x=-10+20*rand(particlesize,narvs); %粒子所在的位置

x=[ 180*rand(particlesize,1)-90 180*2*rand(particlesize,1) 5*rand(particlesize,1)];

v=repmat([3,3,0.3],particlesize,1).*rand(particlesize,3); %粒子的飞翔速度

f=ones(particlesize,1);

FBEST=ones(MaxNum,1); %记录全局最优值得变化

XBEST=ones(MaxNum,narvs); %记录最优值所在位置的变化

for i=1:particlesize

f(i)=POSIfitness(x(i,:));

end

personalbest_x=x;

personalbest_faval=f;

[globalbest_faval,i]=min(personalbest_faval);

globalbest_x=personalbest_x(i,:);

k=1;

while k<=MaxNum

for i=1:particlesize

f(i)=POSIfitness(x(i,:));

if f(i)<personalbest_faval(i) %判断当前位置是否是历史上最佳位置

personalbest_faval(i)=f(i);

personalbest_x(i,:)=x(i,:);

end

[globalbest_faval i]=min(personalbest_faval);

globalbest_x=personalbest_x(i,:);

FBEST(k)=globalbest_faval;

XBEST(k,:)=globalbest_x;

for i=1:particlesize %更新粒子群里每个个体的最新位置

v(i,:)=w*v(i,:)+c1*rand*(personalbest_x(i,:)-x(i,:))...

+c2*rand*(globalbest_x-x(i,:));

for j=1:narvs %判断粒子的飞翔速度是否超过了最大飞翔速度

if v(i,j)>vmax;

v(i,j)=vmax;

elseif v(i,j)<-vmax;

v(i,j)=-vmax;

end

x(i,:)=x(i,:)+v(i,:);

end

if abs(globalbest_faval)<E0

break

end

k=k+1;

end

AGx(q,:)=globalbest_x;

AGfitness(q,:)=globalbest_faval;

end

options=gaoptimset('Generations',100,'Stallgenlimit',10,...

'populationsize',nn,'TolFun',1e-6,'InitialPopulation',AGx, 'plotfcn',@gaplotbestf);

[x,f,reason]=ga(@POSIfitness,narvs,options);

c) Script 3:

function [ f ] = POSIfitness( vars )

%POSIFITNESS 此处显示有关此函数的摘要

% lati 纬度角度制

% longi 经度角度制

lati=vars(1);

longi=vars(2);

L=vars(3);

N=datenum(2015,4,18)-datenum(2015,1,1)+1;

% decli=( 23.45*sin(360*(284+N)/365 /180*pi) )/180*pi;

b=2*pi*N/365;

decli=( 0.006918-0.399912*cos(b)+0.070257*sin(b)-0.006758*cos(2*b)+0.000907*sin(2*b)-0.002697*cos(3*b)+0.00148*sin(3*b) );

%t为北京时间 tt为当地时间即真太阳时

t=(14+42/60): 3/60 :(15+42/60);

tt=t+(longi-120)*4/60;

hour_angle=(tt-12)*15/180*pi;

lati=lati/180*pi;

sin_h=sin(lati)*sin(decli)+cos(lati)*cos(decli)*cos(hour_angle);

% cos_A=( sin(decli)*sin(lati)-cos(decli)*cos(lati)*cos(hour_angle) )./sin_h;

cos_A=( sin_h * sin(lati) - sin(decli) )./ (sqrt(1-sin_h.^2)*cos(lati));

r=L*( sqrt(1-sin_h.^2)./sin_h ); %部分点乘被省略，注意添加

X=[1.0365

1.0699

1.1038

1.1383

1.1732

1.2087

1.2448

1.2815

1.3189

1.3568

1.3955

1.4349

1.4751

1.516

1.5577

1.6003

1.6438

1.6882

1.7337

1.7801

1.8277

];

Y=[0.4973

0.5029

0.5085

0.5142

0.5198

0.5255

0.5311

0.5368

0.5426

0.5483

0.5541

0.5598

0.5657

0.5715

0.5774

0.5833

0.5892

0.5952

0.6013

0.6074

0.6135

];

% y=-r.*cos_A;

% x=r.*sqrt(1-cos_A.^2); %这地方在不同情况下要适时加上-号

% x=r.*cos_A; %最初的方位角理解

% y=r.*sqrt(1-cos_A.^2);

angle1=cos_A(1).*cos_A+sqrt(1-cos_A(1).^2).*sqrt(1-cos_A.^2);

angle2=X(1)/sqrt(X(1)^2+Y(1)^2) * X./sqrt(X.^2+Y.^2) + Y(1)/sqrt(X(1)^2+Y(1)^2) * Y./sqrt(X.^2+Y.^2);

f_a=sum((angle1'-angle2).^2);

f= sum(((X.^2)+(Y.^2)-(r'.^2)).^2 ) + f_a;

% f= sum(((X.^2)+(Y.^2)-(r'.^2)).^2 ) ;

end

d) Script 4:

lati=vars(1);

longi=vars(2);

L=vars(3);

N=datenum(2015,4,18)-datenum(2015,1,1)+1;

% decli=( 23.45*sin(360*(284+N)/365 /180*pi) )/180*pi;

b=2*pi*N/365;

decli=( 0.006918-0.399912*cos(b)+0.070257*sin(b)-0.006758*cos(2*b)+0.000907*sin(2*b)-0.002697*cos(3*b)+0.00148*sin(3*b) );

%t为北京时间 tt为当地时间即真太阳时

t=(14+42/60): 3/60 :(15+42/60);

tt=t+(longi-120)*4/60;

hour_angle=(tt-12)*15/180*pi;

lati=lati/180*pi;

sin_h=sin(lati)*sin(decli)+cos(lati)*cos(decli)*cos(hour_angle);

% cos_A=( sin(decli)*sin(lati)-cos(decli)*cos(lati)*cos(hour_angle) )./sin_h;

cos_A=( sin_h * sin(lati) - sin(decli) )./ (sqrt(1-sin_h.^2)*cos(lati));

r=L*( sqrt(1-sin_h.^2)./sin_h ); %部分点乘被省略，注意添加

X=[1.0365

1.0699

1.1038

1.1383

1.1732

1.2087

1.2448

1.2815

1.3189

1.3568

1.3955

1.4349

1.4751

1.516

1.5577

1.6003

1.6438

1.6882

1.7337

1.7801

1.8277

];

Y=[0.4973

0.5029

0.5085

0.5142

0.5198

0.5255

0.5311

0.5368

0.5426

0.5483

0.5541

0.5598

0.5657

0.5715

0.5774

0.5833

0.5892

0.5952

0.6013

0.6074

0.6135

];

% x=r.*cos_A;

% y=r.*sqrt(1-cos_A.^2);

% f=sum((X-x').^2)+sum((Y-y').^2);

% [x',y']

% angle1=cos_A(1).*cos_A+sqrt(1-cos_A(1).^2).*sqrt(1-cos_A.^2);

% angle2=X(1)/sqrt(X(1)^2+Y(1)^2) * X./sqrt(X.^2+Y.^2) + Y(1)/sqrt(X(1)^2+Y(1)^2) * Y./sqrt(X.^2+Y.^2);

% f_a=sum((angle1-angle2).^2);

% f= sum(((X.^2)+(Y.^2)-(r'.^2)).^2 ) +f_a;

f= sum(((X.^2)+(Y.^2)-(r'.^2)).^2 );

figure

subplot(2,1,1);

plot(t,(X.^2)+(Y.^2))

subplot(2,1,2);

plot(t,r'.^2)

e) Script 5:

fileName = 'F:\学习\课外研究\2015国赛\Appendix4.avi';

obj = VideoReader(fileName);

numFrames = obj.NumberOfFrames;% 帧的总数 61020 每秒5帧

for k = 1 : 2*60*25 : numFrames% 读取数据每两分钟截一次图像共20张

% 对应时间为（8+54/60+6/60/60）:2/60:(9+34/60+46/60/60)

frame = read(obj,k);

imshow(frame);%显示帧

imwrite(frame,strcat(num2str(k),'.jpg'),'jpg');% 保存帧

end

f) Script 6:

rows=ones(21,1);

columes=ones(21,1);

for photos=1:21

fname=[num2str((photos-1)*3) '001.jpg'];

ima=imread(fname);

pic=ima(860:890,1400:1700,:);

thresh=graythresh(pic);

I=im2bw(pic,thresh);

se=strel('diamond',1); % 设置结构元素

temp1=imdilate(I,se);

temp2=imerode(temp1,se);

temp3=noise_ferase(temp2,15,2);

temp4=(logical(temp3-1));

[m,n]=size(temp4);

flags=sum(temp4)-m;

i=1;

tip=flags(n-i+1);

while tip==0

i=i+1;

tip=flags(n-i+1);

end

%倒数过来第i个是影子顶点

colume=n-i+1;

row=floor(mean(find(temp2(:,colume)==0)));

rows(photos)=row;

columes(photos)=colume;

end

%杆长2米占675个像素点以此为参照

Y=(1400-1+columes - 870)/675*2;

X=(860-1+rows - 885)/675*2;

g) Script 7:

function m_new=noise_ferase(b,N,step)

% N 最大的橡皮擦大小，即舍弃上限

% step 橡皮擦每次移动的步长，步长越小，精度越高，但运算速度越慢

%blood=imread('blood001.png');

%w=logical(blood);

[roww,colw]=size(b);

m=b;

m=double(m);m(m==0)=3;m(m==1)=0;m(m==3)=1;

m=logical(m);%黑白颠倒；matlab对应的1为白色、0为黑色，故此时0为边框

m_new=m;

for n=5:5:N %15是边框初值

mm=true(roww+2*n,colw+2*n);

mm(n:n+roww-1,n:n+colw-1)=m_new;

m=mm;

%橡皮擦可以是各种各样的形状，四边形、六边形、八边形等等，这里仅使用四边形，边数越多、越接近于圆，效果越高

%step=2;%探测移动的步长

[rowm,colm]=size(m);

imax=floor((colm-n)/step)+1;

jmax=floor((rowm-n)/step)+1;

for i=1:imax

for j=1:jmax

x=(j-1)*step+1;%首点横坐标

y=(i-1)*step+1;%首点纵坐标

temp1=m(x:x+n-1,y);

temp2=m(x:x+n-1,y+n-1);

temp3=m(x,y:y+n-1)';

temp4=m(x+n-1,y:y+n-1)';

rim=[temp1;temp2;temp3;temp4];%边框之总

flag=all(rim==0);%边框是否包围住噪点

if flag

m(x+1:x+n-2,y+1:y+n-2)=0;

end

m_new=m(n:n+roww-1,n:n+colw-1);

end

imshow(m_new);

h) Script 8:

E0=0.000000001; %允许误差

MaxNum=200; %粒子最大迭代次数

narvs=2; %%%%%%% %目标函数的自变量个数

particlesize=50; %粒子群规模

c1=2; %每个粒子的个体学习因子，也称为加速常数

c2=2; %每个粒子的社会学习因子，也称为加速常数

w=0.75; %惯性因子

vmax=0.8; %粒子的最大飞翔速度

nn=30; % 用于遗传算法的初始种群数量

AGx=ones(nn,narvs);

AGfitness=ones(nn,1);

for q=1:nn

% x=-10+20*rand(particlesize,narvs); %粒子所在的位置

x=[ 180*rand(particlesize,1)-90 180*2*rand(particlesize,1)];

v=repmat([3,3],particlesize,1).*rand(particlesize,narvs); %粒子的飞翔速度

f=ones(particlesize,1);

FBEST=ones(MaxNum,1); %记录全局最优值得变化

XBEST=ones(MaxNum,narvs); %记录最优值所在位置的变化

for i=1:particlesize

f(i)=POSIfitness(x(i,:));

end

personalbest_x=x;

personalbest_faval=f;

[globalbest_faval,i]=min(personalbest_faval);

globalbest_x=personalbest_x(i,:);

k=1;

while k<=MaxNum

for i=1:particlesize

f(i)=POSIfitness(x(i,:));

if f(i)<personalbest_faval(i) %判断当前位置是否是历史上最佳位置

personalbest_faval(i)=f(i);

personalbest_x(i,:)=x(i,:);

end

[globalbest_faval i]=min(personalbest_faval);

globalbest_x=personalbest_x(i,:);

FBEST(k)=globalbest_faval;

XBEST(k,:)=globalbest_x;

for i=1:particlesize %更新粒子群里每个个体的最新位置

v(i,:)=w*v(i,:)+c1*rand*(personalbest_x(i,:)-x(i,:))...

+c2*rand*(globalbest_x-x(i,:));

for j=1:narvs %判断粒子的飞翔速度是否超过了最大飞翔速度

if v(i,j)>vmax;

v(i,j)=vmax;

elseif v(i,j)<-vmax;

v(i,j)=-vmax;

end

x(i,:)=x(i,:)+v(i,:);

end

if abs(globalbest_faval)<E0

break

end

k=k+1;

end

AGx(q,:)=globalbest_x;

AGfitness(q,:)=globalbest_faval;

end

options=gaoptimset('Generations',100,'Stallgenlimit',10,...

'populationsize',nn,'TolFun',1e-6,'InitialPopulation',AGx);%, 'plotfcn',@gaplotbestf);

[x,f,reason]=ga(@POSIfitness,narvs,options);

i) Script 9:

function [ f ] = POSIfitness( vars )

%POSIFITNESS 此处显示有关此函数的摘要

% lati 纬度角度制

% longi 经度角度制

lati=vars(1);

longi=vars(2);

L=2;

N=datenum(2015,7,13)-datenum(2015,1,1)+1;

% decli=( 23.45*sin(360*(284+N)/365 /180*pi) )/180*pi;

b=2*pi*N/365;

decli=( 0.006918-0.399912*cos(b)+0.070257*sin(b)-0.006758*cos(2*b)+0.000907*sin(2*b)-0.002697*cos(3*b)+0.00148*sin(3*b) );

%t为北京时间 tt为当地时间即真太阳时

t=(8+54/60+6/60/60):2/60:(9+34/60+46/60/60);

tt=t+(longi-120)*4/60;

hour_angle=(tt-12)*15/180*pi;

lati=lati/180*pi;

sin_h=sin(lati)*sin(decli)+cos(lati)*cos(decli)*cos(hour_angle);

% cos_A=( sin(decli)*sin(lati)-cos(decli)*cos(lati)*cos(hour_angle) )./sin_h;

cos_A=( sin_h * sin(lati) - sin(decli) )./ (sqrt(1-sin_h.^2)*cos(lati));

r=L*( sqrt(1-sin_h.^2)./sin_h ); %部分点乘被省略，注意添加

X =[ -0.0504

-0.0474

-0.0444

-0.0385

-0.0356

-0.0326

-0.0296

-0.0267

-0.0237

-0.0207

-0.0148

-0.0119

-0.0089

-0.0059

-0.0030

0];

Y =[ 2.3941

2.3644

2.3230

2.2963

2.2726

2.2430

2.2222

2.1896

2.1600

2.1363

2.1067

2.0859

2.0622

2.0267

1.9970

1.9674

1.9407

1.9170

1.8874

1.8696

1.8459];

% y=-r.*cos_A;

% x=r.*sqrt(1-cos_A.^2); %这地方在不同情况下要适时加上-号

% x=r.*cos_A; %最初的方位角理解

% y=r.*sqrt(1-cos_A.^2);

% angle1=cos_A(1).*cos_A+sqrt(1-cos_A(1).^2).*sqrt(1-cos_A.^2);

% angle2=X(1)/sqrt(X(1)^2+Y(1)^2) * X./sqrt(X.^2+Y.^2) + Y(1)/sqrt(X(1)^2+Y(1)^2) * Y./sqrt(X.^2+Y.^2);

% f_a=sum((angle1'-angle2).^2);

% f= sum(((X.^2)+(Y.^2)-(r'.^2)).^2 ) + f_a;

f= sum(((X.^2)+(Y.^2)-(r'.^2)).^2 ) ;

end

III. Internship

As an internship, I was lucky enough that I had a friendly tutor and some intelligent colleagues. They taught me a lot, both in real market and life.

A. Chan theory for trading stocks

1. Introduction:

When I got an internship in GUOTAI JUNAN Future Corp, my daily work was transforming the executive's investment strategies into an algorithm model so that I can use the data to examine the efficacy of the strategy. And my executive, Jianxin Xu, who is a proponent of CHAN theory, insisted that we should contrive a series of rules, by which we can judge when is the right juncture to enter or exit the stock market. But the theory contain a relatively high subjectivity, so we had to adjust some important parameters for different conditions and here is one of those versions.

2. Codes:

a) Script 1:

clear

tic

w=windmatlab;

%load the codes

[w_wss_data,w_wss_codes,w_wss_fields,w_wss_times,w_wss_errorid,w_wss_reqid]...

=w.wss('600000.SH,600004.SH,600005.SH,600006.SH,600007.SH,600008.SH,600009.SH,600010.SH,600011.SH,600012.SH,600015.SH,600016.SH,600017.SH,600018.SH,600019.SH,600020.SH,600021.SH,600022.SH,600023.SH,600026.SH,600027.SH,600028.SH,600029.SH,600030.SH,600031.SH,600033.SH,600035.SH,600036.SH,600037.SH,600038.SH,600039.SH,600048.SH,600050.SH','trade_code');

total_size=size(w_wss_codes,1);

total_balance=nan(total_size,1);

for total_i=1:total_size

codes=w_wss_codes(total_i);

save single_code.mat codes

[w_wsi_data,w_wsi_codes,w_wsi_fields,w_wsi_times,w_wsi_errorid,w_wsi_reqid]...

=w.wsi(codes,'open,high,low,close','2013-07-23 09:00:00','2015-07-22 15:00:00','BarSize=30','Fill=Previous;PriceAdj=F');

data_30ms=[w_wsi_times,w_wsi_data(:,2:3)];

save data_30ms.mat data_30ms

open_prices_30ms=[w_wsi_times,w_wsi_data(:,1)];

save open_prices_30ms.mat open_prices_30ms

[w_wsd_data,w_wsd_codes,w_wsd_fields,w_wsd_times,w_wsd_errorid,w_wsd_reqid]...

=w.wsd(codes,'open,high,low,close','2013-07-23','2015-07-22','Fill=Previous','PriceAdj=F');

data_day=[w_wsd_times,w_wsd_data(:,2:3)];

save data_1day.mat data_day

close_day=[w_wsd_times,w_wsd_data(:,4)];

save data_1day_close.mat close_day

[w_wsd_data,w_wsd_codes,w_wsd_fields,w_wsd_times,w_wsd_errorid,w_wsd_reqid]...

=w.wsd(codes,'open,high,low,close','2013-07-23','2015-07-22','Fill=Previous','PriceAdj=F','Period=W');

data_7days=[w_wsd_times,w_wsd_data(:,2:3)];

save data_7days.mat data_7days

threeK_7days

threeK_1day

threeK_30ms

total_balance(total_i)=balance;

end

save total_balance.mat total_balance

toc

b) Script 2:

% clear

% close

load data_7days.mat

dates=data_7days(:,1); %时间标记

Old=data_7days(:,2:3); %原价格序列【最高价最低价】

N=size(data_7days,1); %z总的序列数量

trend=nan(N,1); %趋势，nan为无，0为反转，1为延续

%Dates 新k线的日期持续范围

%New 新k线的数据

%首先寻找第一个涨跌状态

for i=2:N

flag1=Old(i,1)>Old(i-1,1); %时间i比时间i-1时的最高价的比较

flag2=Old(i,2)>Old(i-1,2); %时间i比时间i-1时的最低价的比较

if flag1==flag2 %出现涨跌状态

n=i; %标记第一个涨跌点

if logical(flag1+flag2) %如果该状态是上升

final_flag=1;

% s_max=max(Old(i,1),Old(i-1,1));

% s_min=max(Old(i,2),Old(i-1,2));

else %如果状态是下降

final_flag=-1;

% s_max=min(Old(i,1),Old(i-1,1));

% s_min=min(Old(i,2),Old(i-1,2));

end

New=Old(i,:);

Dates=[dates(i),dates(i)];

break

end

Enter_weeks=[];

Exit_weeks=[];

for i=n+1:N-1

flag1=Old(i,1)>New(end,1); %时间i比时间i-1时的最高价的比较

flag2=Old(i,2)>New(end,2); %时间i比时间i-1时的最低价的比较

if flag1==flag2 %出现涨跌状态

if logical(flag1+flag2) %如果该状态是上升

new_result=1;

else %如果该状态是下跌

new_result=-1;

end

New=[New;Old(i,:)];

final_flag=[final_flag;new_result];

Dates=[Dates;dates(i),dates(i)];

trend(i)=final_flag(end)==final_flag(end-1);

if final_flag(end-1)==1&&trend(i)==1 %入场信号

Enter_weeks=[Enter_weeks;dates(i+1)]; %入场日期

elseif final_flag(end-1)==-1&&trend(i)==1 %出场信号

Exit_weeks=[Exit_weeks;dates(i+1)]; %出场日期

end

else %未出现涨跌状态，即出现包含状态

if final_flag(end)==-1; %如果上次状态是跌

s_max=min(Old(i,1),Old(i-1,1));

s_min=min(Old(i,2),Old(i-1,2));

else

s_max=max(Old(i,1),Old(i-1,1));

s_min=max(Old(i,2),Old(i-1,2));

end

New(end,:)=[s_max,s_min];

Dates(end,2)=dates(i);

end

temp=size(Enter_weeks,1);

pre_Enter_days=[];

for i=1:temp

pre_Enter_days=[pre_Enter_days;(Enter_weeks(i)-4:Enter_weeks(i))']; %可能成为进场时机的日子

end

temp=size(Exit_weeks,1);

pre_Exit_days=[];

for i=1:temp

pre_Exit_days=[pre_Exit_days;(Exit_weeks(i)-4:Exit_weeks(i))']; %可能成为出场时机的日子

end

c) Script 3:

%clear

%close

load data_1day.mat

dates=data_day(:,1); %时间标记

Old=data_day(:,2:3); %原价格序列【最高价最低价】

N=size(data_day,1); %z总的序列数量

trend=nan(N,1); %趋势，nan为无，0为反转，1为延续

load data_1day_close.mat close_day

n=20;

price_20=nan(N,1); %20天价格

trend_20=nan(N,1); %20天趋势，nan为无，-1为下降，1为上升，0为不变

for i=n-1:N

price_20(i)=mean(close_day(i:i+(n-1),2)); %求出n天均线

end

temp=sign(diff(price_20));

trend_20(n:N)=temp;

%首先寻找第一个涨跌状态

for i=2:N

flag1=Old(i,1)>Old(i-1,1); %时间i比时间i-1时的最高价的比较

flag2=Old(i,2)>Old(i-1,2); %时间i比时间i-1时的最低价的比较

if flag1==flag2 %出现涨跌状态

n=i; %标记第一个涨跌点

if logical(flag1+flag2)&&trend_20(i)==1 %如果该状态是上升.且20天均线是上升

final_flag=1;

% s_max=max(Old(i,1),Old(i-1,1));

% s_min=max(Old(i,2),Old(i-1,2));

else %如果状态是下降

final_flag=-1;

% s_max=min(Old(i,1),Old(i-1,1));

% s_min=min(Old(i,2),Old(i-1,2));

end

New=Old(i,:);

Dates=[dates(i),dates(i)];

break

end

Enter_days=[];

Exit_days=[];

Reduce_days=[];

for i=n+1:N-1

flag1=Old(i,1)>New(end,1); %时间i比时间i-1时的最高价的比较

flag2=Old(i,2)>New(end,2); %时间i比时间i-1时的最低价的比较

if flag1==flag2 %出现涨跌状态

if logical(flag1+flag2)&&trend_20(i)==1 %如果该状态是上升.且20天均线是上升

new_result=1;

else %如果该状态是下跌

new_result=-1;

end

New=[New;Old(i,:)];

final_flag=[final_flag;new_result];

Dates=[Dates;dates(i),dates(i)];

trend(i)=final_flag(end)==final_flag(end-1);

if final_flag(end-1)==1&&trend(i)==1 %入场信号

Enter_days=[Enter_days;dates(i+1)]; %入场日期

end

if final_flag(end-1)==-1&&trend(i)==1 %出场信号止损

Exit_days=[Exit_days;dates(i+1)]; %出场日期

end

if final_flag(end-1)==1&&trend(i)==0 %出场信号止盈

Reduce_days=[Reduce_days;dates(i+1)]; %出场日期

end

else %未出现涨跌状态，即出现包含状态

if final_flag(end)==-1; %如果上次状态是跌

s_max=min(Old(i,1),Old(i-1,1));

s_min=min(Old(i,2),Old(i-1,2));

else

s_max=max(Old(i,1),Old(i-1,1));

s_min=max(Old(i,2),Old(i-1,2));

end

New(end,:)=[s_max,s_min];

Dates(end,2)=dates(i);

end

temp=ismember(Enter_days,pre_Enter_days);

final_Enter_days=Enter_days(temp);

% temp=ismember(Exit_days,pre_Exit_days);

% final_Exit_days=Exit_days(temp);

final_Exit_days=Exit_days; %当趋势为（-1，1）时出场，即该天10：00时的开盘价出场

temp=size(final_Exit_days,1);

temp2=repmat([' 10:00:00'],temp,1);

final_Exit_times=datenum([datestr(final_Exit_days) temp2]);

final_Reduce_days=Reduce_days;

d) Script 4:

%clear

%close

load data_30ms.mat

dates=data_30ms(:,1); %时间标记

Old=data_30ms(:,2:3); %原价格序列【最高价最低价】

load open_prices_30ms.mat %开盘价

open_price=open_prices_30ms(:,2);

N=size(data_30ms,1); %z总的序列数量

trend=nan(N,1); %趋势，nan为无，0为反转，1为延续

load data_1day_close.mat close_day

ddates=close_day(:,1);

dclose=price_20;

temp=size(dates,1);

Temp_p=nan(N,1);

for i=1:temp

temp1= ddates==dates(i);

Temp_p(i)=dclose(temp1); %20天均线

%temp1=find(ddates==Enter_days_for30ms(i));

end

%首先寻找第一个涨跌状态

for i=2:N

flag1=Old(i,1)>Old(i-1,1); %时间i比时间i-1时的最高价的比较

flag2=Old(i,2)>Old(i-1,2); %时间i比时间i-1时的最低价的比较

if flag1==flag2 %出现涨跌状态

n=i; %标记第一个涨跌点

if logical(flag1+flag2)&&Old(i,1)>Temp_p(i) %如果该状态是上升

final_flag=1;

% s_max=max(Old(i,1),Old(i-1,1));

% s_min=max(Old(i,2),Old(i-1,2));

else %如果状态是下降

final_flag=-1;

% s_max=min(Old(i,1),Old(i-1,1));

% s_min=min(Old(i,2),Old(i-1,2));

end

New=Old(i,:);

Dates=[dates(i),dates(i)];

break

end

Enter_times=[];

Reduce_times=[];

for i=n+1:N-1

flag1=Old(i,1)>New(end,1); %时间i比时间i-1时的最高价的比较

flag2=Old(i,2)>New(end,2); %时间i比时间i-1时的最低价的比较

if flag1==flag2 %出现涨跌状态

if logical(flag1+flag2) %如果该状态是上升

new_result=1;

else %如果该状态是下跌

new_result=-1;

end

New=[New;Old(i,:)];

final_flag=[final_flag;new_result];

Dates=[Dates;dates(i),dates(i)];

trend(i)=final_flag(end)==final_flag(end-1);

if final_flag(end-1)==1&&trend(i)==1 %入场信号

Enter_times=[Enter_times;dates(i+1)]; %入场日期

elseif final_flag(end-1)==1&&trend(i)==0 %出场信号止盈

Reduce_times=[Reduce_times;dates(i+1)]; %出场日期

end

else %未出现涨跌状态，即出现包含状态

if final_flag(end)==-1; %如果上次状态是跌

s_max=min(Old(i,1),Old(i-1,1));

s_min=min(Old(i,2),Old(i-1,2));

else

s_max=max(Old(i,1),Old(i-1,1));

s_min=max(Old(i,2),Old(i-1,2));

end

New(end,:)=[s_max,s_min];

Dates(end,2)=dates(i);

end

Enter_days_for30ms=datenum(datestr(Enter_times,'mm/dd/yy')); %进场30分钟线所属于的日子

Reduce_days_for30ms=datenum(datestr(Reduce_times,'mm/dd/yy')); %出场30分钟线所属于的日子

temp1=ismember(Enter_days_for30ms,final_Enter_days);

final_Enter_times=Enter_times(temp1);

temp2=ismember(Reduce_days_for30ms,final_Reduce_days);

final_Reduce_times=Reduce_times(temp2);

Enter_days_for30ms=datenum(datestr(final_Enter_times,'mm/dd/yy')); %合格进场30分钟线所属于的日子

Reduce_days_for30ms=datenum(datestr(final_Reduce_times,'mm/dd/yy')); %合格出场30分钟线所属于的日子

enter_orders=[];

reduce_orders=[];

exit_orders=[];

temp=size(final_Enter_times,1);

for i=1:temp

enter_orders(i)=find(dates==final_Enter_times(i)); %合格的进场时机

end

temp=size(final_Reduce_times,1);

for i=1:temp

reduce_orders(i)=find(dates==final_Reduce_times(i)); %合格的止盈时机

end

temp=size(final_Exit_times,1);

for i=1:temp

try

exit_orders(i)=find(dates==final_Exit_times(i)); %合格的止损时机

catch

continue;

end

exit_orders(exit_orders==0)=[];

Old(enter_orders,:); %进场时的遇到的最好和最低价

Old(reduce_orders,:); %出场时的遇到的最高和最低价

%初始合格点

% plot(dates,open_price);

% hold on

% plot(dates(enter_orders),open_price(enter_orders),'r.')

% plot(dates(exit_orders),open_price(exit_orders),'g.')

% plot(dates(reduce_orders),open_price(reduce_orders),'y.')

%点的筛选

final_orders=nan(N,1);

final_orders(enter_orders)=1; %入场

final_orders(reduce_orders)=0; %减持

final_orders(exit_orders)=-1; %出场

real_enter_orders=[];

real_reduce_orders=[];

real_exit_orders=[];

cash_flow=0;

asset=0;

volume=0;

Account=10^7;

for i=1:N

temp=final_orders(i);

if ~isnan(temp)

if temp==1

if volume<0.50

volume=volume+0.25;

cash_flow=cash_flow-Account*0.25;

asset=asset+Account*0.125/open_price(i);

real_enter_orders=[real_enter_orders;i];

end

elseif temp==0

if volume>0

volume=volume-0.25;

cash_flow=cash_flow+open_price(i)*asset/2;

asset=asset/2;

real_reduce_orders=[real_reduce_orders;i];

end

else

if volume>0

volume=0;

cash_flow=cash_flow+open_price(i)*asset;

asset=0;

real_exit_orders=[real_exit_orders;i];

end

final_close_price=open_prices_30ms(end,2);

balance=asset*final_close_price+cash_flow;

% plot(dates,open_price);

hold off

candle(Old(:,1),Old(:,2),Old(:,1),Old(:,2),[],dates) %烛线图

hold on

uu=8;

plot(dates(real_enter_orders),open_price(real_enter_orders),'r^','markersize',uu,'markerfacecolor','r')

plot(dates(real_exit_orders),open_price(real_exit_orders),'go','markersize',uu,'markerfacecolor','g')

plot(dates(real_reduce_orders),open_price(real_reduce_orders),'ks','markersize',uu,'markerfacecolor','k')

% legend('close_price','buy','exit','reduce','location','best')

xlabel('时间')

ylabel('股价')

load single_code.mat;

names=['股票' cell2mat(codes) '的趋势.png'];

title(['股票' cell2mat(codes) '的趋势']);

saveas(gcf,names)

IV. Projects for my major, Financial Engineering

A. Hedge the risk with options and futures

1. Introduction:

This is the big final project in the instruction Asset Pricing with MATLAB. We were asked to use the options and futures to hedge the risk of any stock portfolio in a proper extent so that we can maximize the utility of customer.

2. Codes:

a) Script 1:

close all

c=yahoo;

%上证50 510050.ss

%上证300 000300.ss

n=2;

codes={'510050.ss','000300.ss'};

index_prices=[];

for i=1:n

stock=fetch(c,codes{i},'close','2014/2/1','2015/2/1');

%雅虎抓取的数据是逆序的，就是最新的数据在最上面

dates=flipud(stock(:,1));

prices=flipud(stock(:,2));

temp1=diff(prices)==0; %价格是否与前一天相同（判断停牌）

temp2=[false;temp1];

prices(temp2)=nan; %停牌当天价格为nan,实际上，因为国内外节假日不同，yahoo把中国

%股市放假时的数据也记录进去了，当然了那时价格也是相对于前一天不变的

%所以我们就也把它当做停牌了对待，删了→。→

fts_temp1=fints(dates,prices);

fts_temp2=fillts(fts_temp1,'linear'); %对停牌数据进行线性插值

fts_temp3=toweekly(fts_temp2); %取周数据,而且这样取得一个好处是，比如国庆那一周没有数据

%但该函数会自动给缺失日期线性插值

prices=fts2mat(fts_temp3);

index_prices=[index_prices,prices];

end

in50=index_prices(:,1);

in300=index_prices(:,2);

n=6;

codes={'600000.ss','600004.ss','600005.ss','600006.ss','600007.ss','600008.ss'};

weight=ones(n,1)/n;

stock=[];

stock_prices=[];

for i=1:n

stock=fetch(c,codes{i},'close','2014/2/1','2015/2/1');

%雅虎抓取的数据是逆序的，就是最新的数据在最上面

dates=flipud(stock(:,1));

prices=flipud(stock(:,2));

temp1=diff(prices)==0; %价格是否与前一天相同（判断停牌）

temp2=[false;temp1];

prices(temp2)=nan; %停牌当天价格为nan,实际上，因为国内外节假日不同，yahoo把中国

%股市放假时的数据也记录进去了，当然了那时价格也是相对于前一天不变的

%所以我们就也把它当做停牌了对待，删了→。→

fts_temp1=fints(dates,prices);

fts_temp2=fillts(fts_temp1,'linear'); %对停牌数据进行线性插值

fts_temp3=toweekly(fts_temp2);%取周数据,而且这样取得一个好处是，比如国庆那一周没有数据

%但该函数会自动给缺失日期线性插值

prices=fts2mat(fts_temp3);

stock_prices=[stock_prices,prices];

end

stockprice=stock_prices*weight; %所购买的股票组合

r_s=diff(log(stockprice)); %股票组合收益率

r_in300=diff(log(in300)); %上证300收益率

r_in50=diff(log(in50)); %上证50收益率

load rf.mat %载入从上海银行间市场网站下载到的无风险利率rf

save hedge_data.mat r_s r_in300 r_in50 rf in50 in300

function [best_hedge_rate,umax]=hedge_option(a)

%% 50ETF期权对冲方法，当前价格S=2.3670，期权的执行价格为2.2，到期期限t=1年，

%输入参数：投资者的风险厌恶系数a

%输出参数：期货的对冲比例n

%因为卖期权需要持有指数，所以再这里我们仅考虑买卖权进行市场风险对冲

%导入数据：

load hedge_data.mat

%计算期权价格所需的参数 r sigma T strike index_price

r=rf(end);

sigma=sqrt(var(r_in50)*52);

%首先求出期权的套期保值比率

T=1; %到期期限

strike=2.2; %执行价格

index_price=in50(end); %指数现值

delta_c=blsdelta(index_price,strike,r,T,sigma); %买权的delta值

delta_p=abs(delta_c-1); %卖权delta值的绝对值

[call_price,put_price]=blsprice(index_price,strike,r,T,sigma); %期权价格

rho=put_price/index_price; %期权价格与指数现价之比，即杠杆率的倒数

%然后对购买的股票组合的超额收益率和市场（上证50）的超额收益率进行回归

C=ones(51,1);

Ri=r_s-rf/52; %股票超额收益率

Rm=r_in50-rf/52; %市场超额收益率

X=[C Rm];

Y=Ri;

[b,bint,r,rint,stats]=regress(Y,X); %回归

beta=b(2); %回归的组合β系数

alpha=b(1); %回归组合的常数项

s=std(r); %残差的标准差——风险体现

N=100;

U=ones(N,1);

nn=linspace(0,beta/delta_p,N); %对冲比例取【0，beta/delta】之间

%当a>>0时，由于投资者是风险厌恶的，所以在beta/delta这点的效用将是最高的，因为这点实现了完全对冲

%若代码换成nn=linspace(0,2*beta/delta_p,N);

%可以看到图形在beta/delta这点处是顶点

%构建股票资产多头和期权空头的投资组合

Port_std=ones(N,1);

for t=1:N

n=nn(t);

PortReturn=((beta-n*delta_p)*mean(Rm)+alpha)/(1+rho*n)+rf(end)/52;%组合的期望收益率

PortRisk=sqrt((beta-n*delta_p)^2*var(Rm)+s^2)/(1+rho*n); %组合的标准差

Port_std(t)=PortRisk;

%效用函数的设定

U(t)=-0.005*a*(PortRisk*100)^2+PortReturn; %a是风险厌恶系数

end

plot(nn,U);

umax=max(U);

best_hedge_rate=nn((U==max(U)));

b) Script 2:

function [best_hedge_rate,umax]=hedge_futures(a)

%% 上证300股指期货对冲方法

%输入参数：投资者的风险厌恶系数a

%输出参数：期货的对冲比例n

%导入数据：

load hedge_data.mat

%首先对股票的超额收益率和市场（上证300）的超额收益率进行回归

C=ones(51,1);

Ri=r_s-rf/52; %股票超额收益率

Rm=r_in300-rf/52; %市场超额收益率

X=[C Rm];

Y=Ri;

[b,bint,r,rint,stats]=regress(Y,X); %回归

beta=b(2); %回归的组合β系数

alpha=b(1); %回归组合的常数项

s=std(r); %残差的标准差——非系统风险的体现

N=100;

U=ones(N,1);

nn=linspace(0,beta,N);%对冲比例取【0，beta】之间

%当a>>0时，由于投资者是风险厌恶的，所以在beta这点的效用将是最高的，因为这点实现了完全对冲

%若代码换成nn=linspace(0,2*beta,N);

%可以看到图形在beta这点处是顶点

%构建股票资产多头和期货空头的投资组合

Port_std=ones(N,1);

for t=1:N

n=nn(t);

PortReturn=((beta-n)*mean(Rm)+alpha)/(1+0.12*n)+rf(end)/52;%组合的期望收益率

PortRisk=sqrt((beta-n)^2*var(Rm)+s^2)/(1+0.12*n);

Port_std(t)=PortRisk;

%效用函数的设定

U(t)=-0.005*a*(PortRisk*100)^2+PortReturn; %a是风险厌恶系数

end

plot(nn,U);

umax=max(U);

best_hedge_rate=nn((U==max(U)));

a=[]; %风险厌恶系数序列

optionhedge=[]; %期权最佳对冲比率

futureshedge=[]; %期货最佳对冲比率

optionutility=[]; %期权最佳对冲下的效用

futuresutility=[]; %期货最佳对冲下的效用

besthedge=[]; %最优对冲策略

bestutility=[]; %最优对冲下的效用

for i=1:101

a(i)=(i-1)*0.005; %风险厌恶系数从0到0.5

[optionhedge(i),optionutility(i)]=hedge_option(a(i));

%特定风险厌恶系数下的期权最优策略

[futureshedge(i),futuresutility(i)]=hedge_futures(a(i));

%特定风险厌恶系数下的期货最优策略

besthedge(i)=max(optionhedge(i),futureshedge(i));

%最优策略的对冲比例

bestutility(i)=max(optionutility(i),futuresutility(i));

%最优策略下的效用

%绘制策略选择图像，红色的表示期权对冲，绿色的表示期货对冲。

if(optionutility(i)>=futuresutility(i))

plot(a(i),optionhedge(i),'ro','markersize',12)

hold on

else

plot(a(i),futureshedge(i),'go','markersize',12)

hold on

end

utility=[futuresutility',optionutility'];

plot(a,utility);

legend('期货对冲策略','期权对冲策略')

title('不同策略的效用图比较')

xlabel('投资者风险厌恶系数')

ylabel('投资效用');

figure

plot(a,besthedge) %最优策略对冲比例图像

figure

plot(a,bestutility) %最优策略下的最大效用

B. Determine parameters in a model

1. Introduction:

Surprisingly, a girl found me through one of my friends and asked a favor. She wanted to figure out some parameters out of complicated equations. It was interesting and I did not adopt the usual way to solve the equations but use the optimum method PSO, particle swarm optimization. I solved it in one day and got a free meal with the girl.

2. Codes:

a) Script 1:

function [ f ] = NSfitness( x )

%x(1) beta0

%x(2) beta1

%x(3) beta2

%x(4) lamda1

%x(5) lamda2

load NSbond.mat;

%载入t r p types

interval=1./types; %付息间隔

firstp=mod(t,interval);%首次付息时间

n=size(p,1);%债券种类个数

p_estimate=ones(n,1);

for i=1:n

if firstp(i)==0

ts=interval(i):interval(i):t(i); %付息时间序列

else

ts=firstp(i):interval(i):t(i); %付息时间序列

end

%%yt=x(1)+x(2)*(1-exp(-x(4)*ts))./(x(4)*ts)+x(3)*((1-exp(-x(5)*ts))./(x(5)*ts)-exp(-x(5)*ts)); %在各付息时间上的利率期限结构

yt=x(1)+x(2)*(1-exp(-x(4)*ts))./(x(4)*ts) + x(3)*((1-exp(-x(4)*ts))./(x(4)*ts)-exp(-x(4)*ts)); %在各付息时间上的利率期限结构

ns=length(ts); %付息期数

cf=ones(ns,1)*(r(i)*interval(i));% 息票

cf(end)=cf(end)+100; %未来现金流

p_estimate(i)=cf'*(exp(-yt.*ts))';

end

%p_err=1/n*sum((p-p_estimate).^2);

p_err=1/n*sum(((p-p_estimate)./dur/sum(dur)).^2); %久期修正

f=p_err;

b) Script 2:

clear all;close all;

tic; %程序运行计时

E0=0.000000001; %允许误差

MaxNum=50; %粒子最大迭代次数

narvs=4; %%%%% %目标函数的自变量个数

particlesize=50; %粒子群规模

c1=2; %每个粒子的个体学习因子，也称为加速常数

c2=2; %每个粒子的社会学习因子，也称为加速常数

w=0.7; %惯性因子

vmax=0.2; %粒子的最大飞翔速度

x=-1+2*rand(particlesize,narvs); %粒子所在的位置

v=0.2*rand(particlesize,narvs); %粒子的飞翔速度

f=ones(particlesize,1);

FBEST=ones(MaxNum,1); %记录全局最优值得变化

XBEST=ones(MaxNum,narvs); %记录最优值所在位置的变化

for i=1:particlesize

f(i)=NSfitness(x(i,:));

end

%f=arrayfun(@NSfitness,x);

personalbest_x=x;

personalbest_faval=f;

[globalbest_faval,i]=min(personalbest_faval);

globalbest_x=personalbest_x(i,:);

k=1;

while k<=MaxNum

for i=1:particlesize

f(i)=NSfitness(x(i,:));

if f(i)<personalbest_faval(i) %判断当前位置是否是历史上最佳位置

personalbest_faval(i)=f(i);

personalbest_x(i,:)=x(i,:);

end

[globalbest_faval i]=min(personalbest_faval);

globalbest_x=personalbest_x(i,:);

FBEST(k)=globalbest_faval;

XBEST(k,:)=globalbest_x;

for i=1:particlesize %更新粒子群里每个个体的最新位置

v(i,:)=w*v(i,:)+c1*rand*(personalbest_x(i,:)-x(i,:))...

+c2*rand*(globalbest_x-x(i,:));

for j=1:narvs %判断粒子的飞翔速度是否超过了最大飞翔速度

if v(i,j)>vmax;

v(i,j)=vmax;

elseif v(i,j)<-vmax;

v(i,j)=-vmax;

end

x(i,:)=x(i,:)+v(i,:);

end

if abs(globalbest_faval)<E0

disp('适应度小于最小误差');

break

end

k=k+1;

end

plot(1:MaxNum,FBEST);

title('全局最优值得变化');

xlabel('粒子群飞翔时间');

legend('算法收敛轨迹')

figure

plot(1:MaxNum,XBEST);

title('最优值位置的变化');

legend('beta0轨迹','beta1轨迹','beta2轨迹','lamda轨迹');

xlabel('粒子群飞翔时间');

toc

save globalbest_x.mat globalbest_x

C. Derived Instruments pricing.

1. Introduction:

In the instruction, Financial Derived Instruments Pricing, professor introduced the principles and methods to price an asset or a derived instrument. So according to the knowledge, I wanted to use different ways to price the asset, including binary tree, differential equation and Monte Carlo simulation.

2. Codes:

a) Script 1:

%%多资产期权风险中性定价——蒙特卡洛模拟

tic

S0 = [50;60;53;73]; %股票初始价格

weight=[0.25;0.25;0.25;0.25]; %每只股票在一篮子资产里所占的权重

stocks_corr=[1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1]; %股票之间的相关系数

Strike=57.5;

n=length(S0);

OptSpec = 'put'; %期权的类型

ValuationDate = '1-Jan-2008'; %期权开始日期

Maturity = '1-Jan-2009'; %期权结束日期

T=(datenum(Maturity)-datenum(ValuationDate))/365; %期权持续期

sigma = [0.20;0.23;0.22;0.30]; %股票年波动率,特别地，这里各只股票之间互不相关

r =0.05; %无风险利率

DividendType = {'continuous'}; %股息形式

q = [0.07;0.05;0.06;0.08]; %股息数量

N=2000;%蒙特卡洛模拟N次

VALUE=ones(N,1); %每次蒙特卡洛模拟得到的解将被放在这里

for i=1:N

M0=mvnrnd(zeros(n,1),stocks_corr); %生成随机数

St=S0.*exp(sigma.*M0.*sqrt(T)+(r-q-sigma.^2/2)*T); %模拟股价运动

St1=S0.*exp(sigma.*(-M0).*sqrt(T)+(r-q-sigma.^2/2)*T);%使用对偶法提高精度

switch OptSpec

case 'call' %看涨期权

Value0=max(weight'*St-Strike,0)*exp(-r*T);

Value1=max(weight'*St1-Strike,0)*exp(-r*T);

Value=0.5*(Value0+Value1);

case 'put' %看跌期权

Value0=max(-weight'*St+Strike,0)*exp(-r*T);

Value1=max(-weight'*St1+Strike,0)*exp(-r*T);

Value=0.5*(Value0+Value1);

end

VALUE(i)=Value;

end

Price=mean(VALUE); %求均值

toc

b) Script 2:

%%方法一：微分方程的公式推导 %买权的公式似乎有点问题。。。

sigma = 0.20; %股票年波动率

r=0.05; %无风险利率

S = 50; %股票初始价格

OptSpec = 'put';

Strike = NaN; %如果Strike是NaN，代表这个期权的执行价格是整个时间段内的最值

%如果Strike是数值，期权的交割价值是整个时间段的最值

ValuationDate = '1-Jan-2008'; %期权开始日期

Maturity = '1-Jan-2009'; %期权结束日期

T=(datenum(Maturity)-datenum(ValuationDate))/365; %期权持续期

q=0.07; %股票的收益率

Fval=[];

switch OptSpec

case 'put' %看跌期权

if isnan(Strike) %执行价格为整个时间段内的最值

M=S; %最值等于初始价格，如果期权已经开始，那么M需要变化！！

a1=(log(M/S)+(-r+q+sigma^2/2)*T)/(sigma*sqrt(T));

a2=a1-(sigma*sqrt(T));

a3=a2+2*(r-q)*sqrt(T)/sigma;

seta=2*(r-q)/sigma^2;

if r~=q

Fval=M*exp(-r*T)*(normcdf(a1)-1/seta*(M/S)^(seta-1)*normcdf(-a3))...

-S*exp(-q*T)*(normcdf(a2)-1/seta*normcdf(-a2));

else

Fval=M*exp(-r*T)*normcdf(a1)...

-S*exp(-q*T)*(normcdf(a2)+log(M/S)*normcdf(-a3)-sigma^2/2*T*normcdf(-a2)...

-sigma*sqrt(T)/sqrt(2*pi)*exp(-a3^2/2));

end

else %最终的交割价值为整个时间段内的最值

disp('该种期权在金融市场上意义很小，并不予以讨论');

end

case 'call' %看涨期权

if isnan(Strike) %执行价格为整个时间段内的最值

M=S; %最值等于初始价格，如果期权已经开始，那么M需要变化！！

b1=(log(M/S)+(r-q+sigma^2/2)*T)/(sigma*sqrt(T));

b2=b1-(sigma*sqrt(T));

b3=-b1+2*(r-q)*sqrt(T)/sigma;

seta=2*(r-q)/sigma^2;

if r~=q

Fval=-M*exp(-r*T)*(normcdf(b2)-1/seta*(M/S)^(seta-1)*normcdf(-b3))...

+S*exp(-q*T)*(normcdf(b1)-1/seta*normcdf(-b1));

else

Fval=-M*exp(-r*T)*normcdf(b2)...

+S*exp(-q*T)*(normcdf(b1)+log(M/S)*normcdf(-b3)-sigma^2/2*T*normcdf(-b1)...

-sigma*sqrt(T)/sqrt(2*pi)*exp(-b3^2/2));

end

else %最终的交割价值为整个时间段内的最值

disp('该种期权在金融市场上意义很小，并不予以讨论');

end

Price=Fval

%%方法二：风险中性定价——蒙特卡洛模拟

%回望期权特征

S0 = 50; %股票初始价格

OptSpec = 'put'; %期权的类型

Strike = NaN; %如果Strike是NaN，代表这个期权的执行价格是整个时间段内的最值

%如果Strike是数值，期权的交割价值是整个时间段的最值

ValuationDate = '1-Jan-2008'; %期权开始日期

Maturity = '1-Jan-2009'; %期权结束日期

Duration=(datenum(Maturity)-datenum(ValuationDate))/365; %期权持续期

Periods=1000; %单位时间个数例如——100时间段 101时间点

t=linspace(0,Duration,Periods+1);

Sigma = 0.20; %股票年波动率

sigma=sqrt(Sigma^2/Periods); %单位时间内的波动率

Rate=0.05; %无风险利率

r=Rate/Periods; %单位时间内的无风险收益率

DividendType = {'continuous'}; %股息形式

DividendAmounts = 0.07; %股息数量

q=DividendAmounts/Periods; %单位时间内的股息率

N=500;%蒙特卡洛模拟N次

VALUE=ones(N,1); %每次蒙特卡洛模拟得到的解将被放在这里

for i=1:N

Mo=randn(Periods+1,1);

MO=cumsum(Mo); %因为收益率服从布朗运动，故不同时间段内的收益率变动相互独立

St=S0*exp(sigma.*MO.*sqrt(t')+(r-q-sigma^2/2).*t'); %模拟股价运动

St1=S0*exp(sigma.*(-MO).*sqrt(t')+(r-q-sigma^2/2).*t');%使用对偶法提高精度

switch OptSpec

case 'call' %看涨期权

if isnan(Strike) %执行价格为整个时间段内的最值

X=min(St);

Value=max(St(end)-X,0)*exp(-Rate*Duration);

X=min(St1);

Value1=max(St1(end)-X,0)*exp(-Rate*Duration);

Value=0.5*(Value+Value1);

else %最终的交割价值为整个时间段内的最值

S=max(St);

Value=max(S-Strike,0)*exp(-Rate*Duration);

S=max(St1);

Value1=max(S-Strike,0)*exp(-Rate*Duration);

Value=0.5*(Value+Value1);

end

case 'put' %看跌期权

if isnan(Strike) %执行价格为整个时间段内的最值

X=max(St);

Value=max(-St(end)*exp(-Rate*Duration)+X,0);

X=max(St1);

Value1=max(-St1(end)*exp(-Rate*Duration)+X,0);

Value=0.5*(Value+Value1);

else %最终的交割价值为整个时间段内的最值

S=min(St);

Value=max(-S+Strike,0)*exp(-Rate*Duration);

S=min(St1);

Value1=max(-S+Strike,0)*exp(-Rate*Duration);

Value=0.5*(Value+Value1);

end

VALUE(i)=Value;

end

Price=mean(VALUE) %求均值

%%方法三：二叉树定价

%根据衍生金融工具箱

%回望期权特征

OptSpec = 'put'; %看涨还是看跌

Strike = NaN; %如果Strike是NaN，代表这个期权的执行价格是整个时间段内的最值

%如果Strike是数值，期权的交割价值是整个时间段的最值

Settle = '01-Jan-2008'; %期权开始日期

ExerciseDates = '01-Jan-2009'; %期权结束日期

InstSet = instlookback(OptSpec, Strike, Settle, ExerciseDates); %设置期权特征结构

%股票特征

AssetPrice = 50; %股票现价

Sigma = 0.20; %股票波动性

DividendType = {'continuous'}; %股票分红形式

DividendAmounts = 0.07; %股票红利率

StockSpec = stockspec(Sigma, AssetPrice, DividendType, DividendAmounts);%, ExDividendDates);

%设置股票结构特征

%时间特征

Periods=100; %二叉树的期数，期数越大，得出的结果越精确

ValuationDate = {'1-Jan-2008'}; %开始时间

Maturity = '1-Jan-2009'; %结束时间

TimeSpec = crrtimespec(ValuationDate, Maturity, Periods); %设置时间结构特征

%利率特征

StartDates = '01-Jan-2008';

EndDates = '01-Jan-2009';

Rates=0.05; %无风险利率

RateSpec = intenvset('Rates', Rates, 'StartDates',StartDates,...

'EndDates', EndDates, 'Compounding', 365); %设置无风险利率期限结构特征

%股票二叉树

CRRTree = crrtree(StockSpec, RateSpec, TimeSpec); %由上得出股票二叉树

%定价

[Price, PriceTree] = crrprice(CRRTree, InstSet); %由股票二叉树和期权特征得出定价

%或者 Price = lookbackbycrr(CRRTree, OptSpec, Strike, Settle,ExerciseDates);

Price

c) Script 3:

%3年期固定息票支付，息票率7%，本金100，每年支付2次，天数计算规则Basic=1

%利率期限结构

%%时间结构 hjmtimespec

Compounding=2;

ValuationDate='01-01-2010';

StartDates='01-01-2010';

Maturity=datemnth('01-01-2010',6*(1:6)');

TimeSpec=hjmtimespec(ValuationDate,Maturity,Compounding);

%%利率结构 hjmratespec

Rates=[0.04,0.05,0.058,0.064,0.069,0.073]';

Basis=1;

Maturity=datemnth('01-01-2010',5*(1:6)');

RateSpec=intenvset('Rates',Rates,'Compounding',Compounding,...

'StartDates',StartDates,'EndDates',Maturity,'Basis',Basis);

%%利率波动结构 hjmvolspec

%Volatility = [.02; .019; .018; .017; .016; .018];

%CurveTerm = [1; 2; 3; 4; 5 ;6];

%VolSpec = hjmvolspec('Stationary', Volatility , CurveTerm);

Sigma_0=0.01;

VolSpec = hjmvolspec('Constant', Sigma_0);

%%利率树

HJMTree=hjmtree(VolSpec,RateSpec,TimeSpec);

treeviewer(HJMTree);

%%产品特征 hjminstset

CouponRate=0.07;

Settle=datenum('01-01-2010');

Reset=2;

Principal=100;

EndMonthRule=1;

InstSet=instfixed(CouponRate, Settle, Maturity(end), Reset, Basis, Principal,EndMonthRule);

%%汇总定价

[Price,PriceTree]=hjmprice(HJMTree,InstSet);

treeviewer(PriceTree)

D. CDOs Pricing

1. Introduction:

This is a project in Financial Risk Management. Here we were going to analyze the risk of CMO, a derived financial instrument basing on the elementary mortgage loan. Hence we designed a model and regarded the interest rate as the pivot ingredient to calculate the forthcoming cash current. In the meanwhile we used CVaR to estimate prepayment risk.

2. Codes:

a) Script 1:

%r=0:0.001:0.08

%y=arrayfun(@FRM_CMO,arrayfun(@cal_PSA,r))

clear

tic

PAC=92167402;

A=212846852;

C=143576072;

Z=131292782;

F0=[PAC,A,C,Z];

n=100;

dur=ones(n,4); %各现金流久期

value=ones(n,4); %四个现金流的折现价值，折现率使用的是无风险利率对应的利率期限结构

%个人认为对冲基金会把风险对冲完整，故而对应的是风险中性的世界

para=linspace(-0.0277,0.0277,n); %利率的变动，0.0227是算出来的利率期限结构中一个月期的利率强度

for i=1:100

[value(i,:),dur(i,:)]=FRM_CMO(para(i));

end

for i=1:4

figure

plot(0.0277+para,value(:,i));

title(['第' num2str(i) '档产品随价格随利率的变动']);

xlabel('利率');

ylabel('价值');

axis tight

end

n=1000;

dur=ones(n,4); %各现金流久期

total_value=ones(n,4); %四个现金流的加权久期

value=ones(n,4); %四个现金流的折现价值，折现率使用的是无风险利率对应的利率期限结构

%个人认为对冲基金会把风险对冲完整，故而对应的是风险中性的世界

para=normrnd(0,0.005,n,1); %利率的变动

for i=1:n

[value(i,:),dur(i,:),total_value(i,:)]=FRM_CMO(para(i));

end

CMO_VAR; %计算VAR值或者CMO_VAR2

cf_rmean=(mean(value)-F0);%./F0; %平均收益率

cf_std=std(value);%./mean(value); %标准化后的波动率

cf_var=cfVAR;%./F0; %var相对于本金的占比，视作另一种风险度量

CV=cf_std./cf_rmean;

CVvar=cf_var./cf_rmean;

[temp1,temp2,temp3]=FRM_CMO(0);

RAROC=(temp3)./cf_var;

value_0=FRM_CMO(0);

value_1=FRM_CMO(0.005);

value_2=FRM_CMO(-0.005);

eDur=-(value_1-value_2)/0.01./value_0; %有效久期=（利率上升0.5%时的价值-利率下降0.5%时的价值）/1%/现值

toc

b) Script 2:

function [ABCD_value,Dur,total_value]=FRM_CMO(para)

%para为利率的浮动大小

F0=579883108;

Terms=357;

WAC=0.0675;

PT=0.0520;

PAC=92167402;

A=212846852;

C=143576072;

Z=131292782;

%加入利率期限结构

ts=1/12:1/12:Terms/12;

r=[2.890,3.1199,3.5728,3.6973,3.9190,4.2665,4.5958];

r1=r/100; %原利率结构

r2=r/100+para; %新利率结构

t=[3/12,6/12,2,3,5,10,30];

yield=interp1(t,r,ts,'PCHIP')/100;

yield=yield+para; %插值

Dur=ones(1,4);

discount=(1+yield).^(-ts);

discount1=(1+r1).^(-t);

discount2=(1+r2).^(-t);

load CMO_weight.mat %载入零息债券的权重

bondvalue=[PAC*x(:,1),A*x(:,2),C*x(:,3),Z*x(:,4)];

Z_Bond=(1./discount1.*discount2)*bondvalue;

PSA=cal_PSA(yield(1));

%PSA=1.98; %缺少利率期限结构的输入

CPR=min([0.002:0.002:Terms*0.002]',0.06);

SMM=1-(1-PSA*CPR).^(1/12);

ts=(1:Terms)';

F=ones(Terms,1); %月初本金余额循环到最后，这个会多出一个来

F(1)=F0;

payment=ones(Terms,1); %月偿付额

I=ones(Terms,1); %MBS的利息

PI=ones(Terms,1); %CMO的利息

pay_sche=ones(Terms,1); %计划偿还本金

pay_pre=ones(Terms,1); %提前偿付本金

CF=ones(Terms,1); %流向CMO的总现金流

for t=1:Terms

payment(t)=F(t)/((1-1/(1+WAC/12)^(Terms-t+1))/(WAC/12));

I(t)=F(t)*(WAC/12);

PI(t)=F(t)*(PT/12);

pay_sche(t)=payment(t)-I(t);

pay_pre(t)=F(t)*SMM(t);

CF(t)=PI(t)+pay_sche(t)+pay_pre(t);

F(t+1)=F(t)-(pay_sche(t)+pay_pre(t));

end

pay_prin=pay_sche+pay_pre; %CMO总本金偿付

F(end)=[];

Amat=[ts,F,payment,PI,I,pay_sche,pay_pre,CF];

AAmat=['时间','月初本金','月偿还额','MBS利息','CMO利息','计划本金偿还','提前本金偿还','CMO现金流'];

load PAC_p.mat %载入PAC_ppay

PAC_prin=ones(Terms,1)*PAC; %期初的本金

A_prin=ones(Terms,1)*A;

C_prin=ones(Terms,1)*C;

Z_prin=ones(Terms,1)*Z;

PAC_ipay=zeros(Terms,1);

A_ppay=zeros(Terms,1);

A_ipay=zeros(Terms,1);

C_ppay=zeros(Terms,1);

C_ipay=zeros(Terms,1);

Z_ppay=zeros(Terms,1);

Z_ipay=zeros(Terms,1);

cf=CF;

%CF是每一期的现金流总支付

for i=1:Terms

PAC_ipay(i)=PAC_prin(i)*PT/12;

A_ipay(i)=A_prin(i)*PT/12;

C_ipay(i)=C_prin(i)*PT/12;

Z_ipay(i)=Z_prin(i)*PT/12; %%%

Z_prin(i+1)=Z_prin(i)+Z_ipay(i); %%% %Z 不接收利息

Z_ipay(i)=0;

prin_pay=cf(i)-PAC_ipay(i)-A_ipay(i)-C_ipay(i); %扣除的部分不含Z的利息

PAC_prin(i+1)=PAC_prin(i)-PAC_ppay(i);

prin_pay=prin_pay-PAC_ppay(i);

A_ppay(i)=prin_pay;

A_prin(i+1)=A_prin(i)-prin_pay;

if A_prin(i+1)<=0

A_prin(i+1)=0;

A_ppay(i)=A_prin(i);

prin_pay=prin_pay-A_prin(i);

C_ppay(i)=prin_pay;

C_prin(i+1)=C_prin(i)-prin_pay;

if C_prin(i+1)<=0

C_prin(i+1)=0;

C_ppay(i)=C_prin(i);

prin_pay=prin_pay-C_prin(i);

Z_ppay(i)=prin_pay;

Z_prin(i+1)=Z_prin(i+1)-prin_pay;

end

%[PAC_ppay,A_ppay,C_ppay,Z_ppay]

ABCD_prin=[PAC_prin,A_prin,C_prin,Z_prin];

ABCD_ppay=[PAC_ppay,A_ppay,C_ppay,Z_ppay];

ABCD_ipay=[PAC_ipay,A_ipay,C_ipay,Z_ipay];

ABCD_cf=ABCD_ppay+ABCD_ipay;

ts=1/12:1/12:Terms/12;

ABCD_value=ones(1,4);

for i=1:4

ABCD_value(i)=(ABCD_cf(:,i)'*discount'); %%%在风险中性世界中的价值

Dur(i)=(ts'.*ABCD_cf(:,i))'*discount'/(ABCD_cf(:,i)'*discount');

end

total_value=ABCD_value - Z_Bond;

totalDur=Dur*([PAC;A;C;Z])/sum([PAC;A;C;Z]);

%Dur;

%r=0:0.001:0.08

%y=arrayfun(@FRM_CMO,arrayfun(@cal_PSA,r))

c) Script 3:

%r=0:0.001:0.08

%y=arrayfun(@FRM_CMO,arrayfun(@cal_PSA,r))

clear

tic

PAC=92167402;

A=212846852;

C=143576072;

Z=131292782;

F0=[PAC,A,C,Z];

n=100;

dur=ones(n,4); %各现金流久期

value=ones(n,4); %四个现金流的折现价值，折现率使用的是无风险利率对应的利率期限结构

%个人认为对冲基金会把风险对冲完整，故而对应的是风险中性的世界

para=linspace(-0.0277,0.0277,n); %利率的变动，0.0227是算出来的利率期限结构中一个月期的利率强度

for i=1:100

[value(i,:),dur(i,:)]=FRM_CMO(para(i));

end

for i=1:4

figure

plot(0.0277+para,value(:,i));

title(['第' num2str(i) '档产品随价格随利率的变动']);

xlabel('利率');

ylabel('价值');

axis tight

end

n=1000;

dur=ones(n,4); %各现金流久期

total_value=ones(n,4); %四个现金流的加权久期

value=ones(n,4); %四个现金流的折现价值，折现率使用的是无风险利率对应的利率期限结构

%个人认为对冲基金会把风险对冲完整，故而对应的是风险中性的世界

para=normrnd(0,0.005,n,1); %利率的变动

for i=1:n

[value(i,:),dur(i,:),total_value(i,:)]=FRM_CMO(para(i));

end

CMO_VAR; %计算VAR值或者CMO_VAR2

cf_rmean=(mean(value)-F0);%./F0; %平均收益率

cf_std=std(value);%./mean(value); %标准化后的波动率

cf_var=cfVAR;%./F0; %var相对于本金的占比，视作另一种风险度量

CV=cf_std./cf_rmean;

CVvar=cf_var./cf_rmean;

[temp1,temp2,temp3]=FRM_CMO(0);

RAROC=(temp3)./cf_var;

value_0=FRM_CMO(0);

value_1=FRM_CMO(0.005);

value_2=FRM_CMO(-0.005);

eDur=-(value_1-value_2)/0.01./value_0; %有效久期=（利率上升0.5%时的价值-利率下降0.5%时的价值）/1%/现值

toc

E. Efficient Frontier

1. Introduction:

This is the classic model CAPM. Everybody want to get returns with the lowest volatility. And here is a way to atain the goal under model CAPM.

2. Codes:

a) Script 1:

close all

c=yahoo;

n=6;

codes={'600629.ss','000402.sz','600383.ss','600018.ss','600361.ss','600785.ss'};

stock=[];

stock_prices=[];

newopen=ones(1,n);

for i=1:n

stock=fetch(c,codes{i},{'close','open'},'2014/4/13','2015/6/13');

%雅虎抓取的数据是逆序的，就是最新的数据在最上面

%newopen(i)=stock(1,3);

dates=flipud(stock(:,1));

prices=flipud(stock(:,2));

temp1=diff(prices)==0; %价格是否与前一天相同（判断停牌）

temp2=[false;temp1];

prices(temp2)=nan; %停牌当天价格为nan,实际上，因为国内外节假日不同，yahoo把中国

%股市放假时的数据也记录进去了，当然了那时价格也是相对于前一天不变的

%所以我们就也把它当做停牌了对待，删了→。→

stock_prices=[stock_prices,prices];

end

%如果连不上，就直接用 load stockdata.mat

ExpReturn=nanmean(price2ret(stock_prices)); %nanmean剔除nan数据后的估计

ExpCovariance=nancov(price2ret(stock_prices));

NumPorts=100;%资产组合有效前沿上取的100个点

BorrowRate=0.10/365;

RisklessRate=0.02/365; %假定无风险利率为0.02，那么转变成日数据就是0.02/365

RiskAversion=3;

b) Script 2:

%% 有卖空限制

[PortRisk,PortReturn,PortWts]=portopt(ExpReturn,ExpCovariance,NumPorts);%有效前沿上的点

plot(PortRisk,PortReturn,'r');

plot(PortRisk(1),PortReturn(1),'ro') %最小方差点

cv=(PortReturn-BorrowRate)./PortRisk;

index=find(cv==max(cv));

plot(PortRisk(index),PortReturn(index),'r*') %最佳风险资产组合点

p=[max(cv),BorrowRate];

CML=polyval(p,PortRisk);

plot(PortRisk,CML,'r--') %CML线

[RiskyRisk,RiskyReturn,RiskyWts,RiskyFraction,OverallRisk,OverallReturn]=portalloc(PortRisk,PortReturn,PortWts,RisklessRate,BorrowRate,RiskAversion)

plot(RiskyRisk,RiskyReturn,'r^') %最佳组合点

legend('无卖空限制条件下的有效前沿','最小方差点','最佳风险组合','CML线','最佳组合点',...

'有卖空限制条件下的有效前沿','最小方差点','最佳风险组合','CML线','最佳组合点','location','best')

V. Projects applied in stock market

A. My idea for picking the most profitable stock

1. Introduction:

The frantic stock market in the early 2015 made people exciting, all of them wanted to make a fortune in this trend. So did I. And I believed that the financial statements of companies would tell the fact whether the companies were healthy enough, eventually we could get the conclusion whether it was worthwhile to buy a specific company. But the question is, what could be the standard for a healthy and potential company. So here I use the classic method PCA, principle component analysis, to calculate the best company.

2. Codes:

% PCA法选择最优企业

load PCAdata.mat %处理好的数据，均是效益型指标

Data;

std_indic=zscore(Data);

R=corrcoef(std_indic);

[vec1,lamda,rate]=pcacov(R);

f=repmat(sign(sum(vec1)),size(vec1,1),1);

vec2=vec1.*f;

per=0.85;

num=find(cumsum(rate/100)>per,1);

p_score=std_indic*vec2(:,1:num);

total_score=p_score*rate(1:num);

[sort_ts,ind]=sort(total_score,'descend'); % 第ind(i)只是最好的企业

B. The incongruity between high return and high risk

1. Introduction:

China stock market got a crazy boost in the early 2015. But Professor.Zhang Xiaojun told us something interesting about a 1990s paper, whose theme is that the low volatility will bring a long-term high return. And I was curious about its authenticity so that I programmed and took the Chinese data into the test.

2. Codes:

%股票数据下载

close all

clear

load StockCodes.mat

tic

% n=size(Codes,1);

n=1032;

% a=ones(n,1)

% mat2cell(Codes,a)

flags=ones(n,1);%指示抓取数据是否成功

s=struct();

for i=1:n

try

c=yahoo;

stock=fetch(c,Codes(i,:),'close','2013/1/1','2015/1/1');

%雅虎抓取的数据是逆序的，就是最新的数据在最上面

dates=flipud(stock(:,1));

prices=flipud(stock(:,2));

temp1=diff(prices)==0; %价格是否与前一天相同（判断停牌）

temp2=[false;temp1];

prices(temp2)=nan; %停牌当天价格为nan,实际上，因为国内外节假日不同，yahoo把中国

%股市放假时的数据也记录进去了，当然了那时价格也是相对于前一天不变的

%所以我们就也把它当做停牌了对待，删了→。→

fts_temp1=fints(dates,prices);

fts_temp2=fillts(fts_temp1,'linear'); %对停牌数据进行线性插值

fts_temp3=toweekly(fts_temp2);%取周数据,而且这样取得一个好处是，比如国庆那一周没有数据

%但该函数会自动给缺失日期线性插值

prices=fts2mat(fts_temp3);

s(i).prices=prices;

s(i).dates=fts_temp3.dates;

disp(['抓取第 ' num2str(i) ' 只股票成功'])

catch %抓取失败，进行标记

flags(i)=0;

disp(['抓取第 ' num2str(i) ' 只股票失败'])

end

% i

end

toc %抓到852只股票，用时2900s

load sss.mat

nums=sum(flags);

orders=find(flags);

counters=0;

load rf_3m.mat

rf=rf_3m/100/52; %转化成周数据

for i=1:nums

j=orders(i);

if size(s(j).prices,1)<90|sum(isnan(s(j).prices))>10

disp(['第 ' num2str(j) ' 只股票数据太少']);

counters=counters+1;

else

r_s=diff(log(s(j).prices)); %股票周收益率

r_s_total(i)=nansum(r_s); %股票这两年来的总收益率

covr=nancov(r_s-rf,r_index-rf);

beta(i)=covr(1,2)/nanvar(r_index-rf);

v_s(i)=nanstd((r_s)); %股票周波动率

end

fs=[v_s',r_s_total',beta'];

%16只股票数据太少，剩余852只

%总风险

temp1=isnan(fs); %去nan

temp2=logical(sum(temp1,2));

fs(temp2,:)=[];

temp1=(fs==0); %去0

temp2=logical(sum(temp1,2));

fs(temp2,:)=[];

fs1=sortrows(fs);

%834

nn=40;

plot(fs1(1:nn,1),fs1(1:nn,2),'b.');

mean(fs1(1:nn,2))

hold on

plot(fs1(end-nn+1:end,1),fs1(end-nn+1:end,2),'r.');

mean(fs1(end-nn+1:end,2))

hold off

legend('低风险股票','高风险股票')

xlabel('股票周波动率')

ylabel('股票长期收益')

title('股票长期收益与其波动性的关系')

%与系统风险

figure

fs2=[(fs(:,3)),fs(:,2)];

temp=fs2(:,1)>0;

fs2=fs2(temp,:);

[fs1,I]=sortrows(fs2);

%834

% beta=fs(:,3);

% beta=beta(I);

nn=80;

plot(fs1(1:nn,1),fs1(1:nn,2),'b.');

mean(fs1(1:nn,2))

hold on

plot(fs1(end-nn+1:end,1),fs1(end-nn+1:end,2),'r.');

mean(fs1(end-nn+1:end,2))

hold off

legend('低beta值股票','高beta值股票')

xlabel('股票周波动率')

ylabel('股票长期收益')

title('股票长期收益与其系统风险的关系')

figure

plot(fs1(:,1),fs1(:,2),'b.');

C. Granville's thoughts

1. Introduction:

Granville's eight rules in the stock market are the classic stock-picking methods in technical analysis. When professor taught us how to trade stocks in the class, he admonished that we should attach importance to the price moving curve as Granville once said. As a fresh man in the stock market, however, I did not believe these casual curve would include some crucial imformations. I doubted the theory so I designed the codes to examine it.

2. Codes:

a) Script 1:

function flag=GLB_prin(StockCode,StartDate,EndDate)

try

c=yahoo;

b=fetch(c,StockCode,'close',StartDate,EndDate); %抓取数据

n=20; %n天移动平均

psize=length(b); %数据时间长度

b=flipud(b); %将数据倒回来

dates=b(:,1); %时间

bmean=ones(psize-(n-1),1); %n天均线

bfints=fints(b);

for i=1:psize-(n-1)

bmean(i)=mean(b(i:i+(n-1),2)); %求出n天均线

end

bb=b(n:end,:); %日线

bbmean=[dates(n:end),bmean]; %n天均线

bmfints=fints(bbmean);

%以上画图显示股价和20天均线

newdates=dates(n:end);

impert=(bb(:,2)-bbmean(:,2))./bbmean(:,2);%乖离率

%法则一，今天和前一天进行比较，收盘价

diffbmean=diff(bmean); %判断20天均线是否上升

meanflag=diffbmean>0; %移动平均线相较于前一天是上升时，meanflag为1，反之为0；

bmflag=bbmean(:,2)>bb(:,2);%20天平均线是否高于股价，若高于，则bmflag为1，反之为0；

diffbm=diff(bmflag);%若20天平均线由低于股价变为高于股价，则diffbm为1，即20天移动平均线上穿股价

%反之-1，代表20天移动平均线下穿股价

GLB1flag=diffbm==-1&meanflag==1; %判断信号

GLB1flag=logical([0;GLB1flag]); %补齐数据

GLB1posi=newdates(GLB1flag); %信号时间

GLB1price=bb(GLB1flag,2); %信号价格

GLB1fints=fints([GLB1posi,GLB1price]);

%我们认为2法则和1法则一样，故跳过

%法则三,今天和前一天进行比较，在今天决策

GLB3_1=impert<0.02&impert>=0; GLB3_1(end)=[];

GLB3_2=meanflag==1;

GLB3_3=diff(impert)>0.02;

GLB3flag=GLB3_1&GLB3_2&GLB3_3;

GLB3flag=logical([0;GLB3flag]);

GLB3posi=newdates(GLB3flag);

GLB3price=bb(GLB3flag,2);

GLB3fints=fints([GLB3posi,GLB3price]);

%法则四，乖离率一过大就买入

GLB4_1=impert<-0.07; GLB4_1(end)=[];

GLB4_2=meanflag==0;

%GLB4_3=diff(impert)>0.02;

GLB4flag=GLB4_1&GLB4_2;

GLB4flag=logical([0;GLB4flag]);

GLB4posi=newdates(GLB4flag);

GLB4price=bb(GLB4flag,2);

GLB4fints=fints([GLB4posi,GLB4price]);

%法则五，今天和前一天进行比较，收盘价

GLB5flag=diffbm==1&meanflag==0;

GLB5flag=logical([0;GLB5flag]);

GLB5posi=newdates(GLB5flag);

GLB5price=bb(GLB5flag,2);

GLB5fints=fints([GLB5posi,GLB5price]);

%我们认为6法则和5法则一样，故跳过

%法则七,今天和前一天进行比较，在今天决策

GLB7_1=impert>-0.02&impert<=0; GLB7_1(end)=[];

GLB7_2=meanflag==0;

GLB7_3=diff(impert)<-0.02;

GLB7flag=GLB7_1&GLB7_2&GLB7_3;

GLB7flag=logical([0;GLB7flag]);

GLB7posi=newdates(GLB7flag);

GLB7price=bb(GLB7flag,2);

GLB7fints=fints([GLB7posi,GLB7price]);

%法则八，乖离率一过大就买入

GLB8_1=impert>0.07; GLB8_1(end)=[];

GLB8_2=meanflag==1;

%GLB4_3=diff(impert)>0.02;

GLB8flag=GLB8_1&GLB8_2;

GLB8flag=logical([0;GLB8flag]);

GLB8posi=newdates(GLB8flag);

GLB8price=bb(GLB8flag,2);

GLB8fints=fints([GLB8posi,GLB8price]);

%计算盈亏情况

flag1=GLB1flag*1+GLB3flag*1+GLB4flag*1-(GLB5flag*1+GLB7flag*1+GLB8flag*1); %决定买卖点

flag2=sign(flag1);

temp1=diff(flag2)==0; %找出无效的连续买点信号和卖点信号

temp2=[false;temp1]; %因为差分消除了一个数据，故而此处补上一个0

flag2(temp2)=0; %消除无效信号

q=flag2((flag2~=0)); %提取出带有信号时间点

qix=find(flag2~=0); %标记这些点的位置

temp1=diff(q)==0; %如果这些时间点是连续同方向操作，视为无效信号

temp2=[false;temp1]; %补齐差分所消数据

q(temp2)=0; %消除无效信号

finalflag=zeros(size(newdates)); %1为买进信号，-1为卖出信号

finalflag(qix(q==1))=1; %在时间序列上标记最后买卖点信号

finalflag(qix(q==-1))=-1;

buytime=qix(q==1); %买点时间

soldtime=qix(q==-1); %卖点时间

flag=finalflag(end); %最后一天是否是买点、卖点

catch

flag=nan;

end

b) Script 2:

%b浦发银行

close;

clear;

load stockdata.mat %fetch寻找的数据

n=20; %n天移动平均

psize=length(b);

b=flipud(b);

dates=b(:,1);

bmean=ones(psize-(n-1),1);

bfints=fints(b);

plot(bfints,'k');

for i=1:psize-(n-1)

bmean(i)=mean(b(i:i+(n-1),2));

end

hold on

bb=b(n:end,:);

bbmean=[dates(n:end),bmean];

bmfints=fints(bbmean);

plot(bmfints,'r');

%以上画图显示股价和20天均线

newdates=dates(n:end);

impert=(bb(:,2)-bbmean(:,2))./bbmean(:,2);%乖离率

%法则一，今天和前一天进行比较，收盘价

diffbmean=diff(bmean);

meanflag=diffbmean>0;%移动平均线相较于前一天是上升时，meanflag为1，反之为0；

bmflag=bbmean(:,2)>bb(:,2);%20天平均线是否高于股价，若高于，则bmflag为1，反之为0；

diffbm=diff(bmflag);%若20天平均线由低于股价变为高于股价，则diffbm为1，即20天移动平均线上穿股价

%反之-1，代表20天移动平均线下穿股价

GLB1flag=diffbm==-1&meanflag==1;

GLB1flag=logical([0;GLB1flag]);

GLB1posi=newdates(GLB1flag);

GLB1price=bb(GLB1flag,2);

GLB1fints=fints([GLB1posi,GLB1price]);

plot(GLB1fints,'b.');%葛兰碧一法则

%我们认为2法则和1法则一样，故跳过

%法则三,今天和前一天进行比较，在今天决策

GLB3_1=impert<0.02&impert>=0; GLB3_1(end)=[];

GLB3_2=meanflag==1;

GLB3_3=diff(impert)>0.02;

GLB3flag=GLB3_1&GLB3_2&GLB3_3;

GLB3flag=logical([0;GLB3flag]);

GLB3posi=newdates(GLB3flag);

GLB3price=bb(GLB3flag,2);

GLB3fints=fints([GLB3posi,GLB3price]);

plot(GLB3fints,'m.');

%法则四，乖离率一过大就买入

GLB4_1=impert<-0.07; GLB4_1(end)=[];

GLB4_2=meanflag==0;

%GLB4_3=diff(impert)>0.02;

GLB4flag=GLB4_1&GLB4_2;

GLB4flag=logical([0;GLB4flag]);

GLB4posi=newdates(GLB4flag);

GLB4price=bb(GLB4flag,2);

GLB4fints=fints([GLB4posi,GLB4price]);

plot(GLB4fints,'g.');

%法则五，今天和前一天进行比较，收盘价

GLB5flag=diffbm==1&meanflag==0;

GLB5flag=logical([0;GLB5flag]);

GLB5posi=newdates(GLB5flag);

GLB5price=bb(GLB5flag,2);

GLB5fints=fints([GLB5posi,GLB5price]);

plot(GLB5fints,'bo');%葛兰碧一法则

%我们认为6法则和5法则一样，故跳过

%法则七,今天和前一天进行比较，在今天决策

GLB7_1=impert>-0.02&impert<=0; GLB7_1(end)=[];

GLB7_2=meanflag==0;

GLB7_3=diff(impert)<-0.02;

GLB7flag=GLB7_1&GLB7_2&GLB7_3;

GLB7flag=logical([0;GLB7flag]);

GLB7posi=newdates(GLB7flag);

GLB7price=bb(GLB7flag,2);

GLB7fints=fints([GLB7posi,GLB7price]);

plot(GLB7fints,'mo');

%法则八，乖离率一过大就买入

GLB8_1=impert>0.07; GLB8_1(end)=[];

GLB8_2=meanflag==1;

%GLB4_3=diff(impert)>0.02;

GLB8flag=GLB8_1&GLB8_2;

GLB8flag=logical([0;GLB8flag]);

GLB8posi=newdates(GLB8flag);

GLB8price=bb(GLB8flag,2);

GLB8fints=fints([GLB8posi,GLB8price]);

plot(GLB8fints,'go');

%计算盈亏情况

flag1=GLB1flag*1+GLB3flag*1+GLB4flag*1-(GLB5flag*1+GLB7flag*1+GLB8flag*1);

flag2=sign(flag1);

temp1=diff(flag2)==0;

temp2=[false;temp1];

flag2(temp2)=0;

q=flag2((flag2~=0));

qix=find(flag2~=0);

temp1=diff(q)==0;

temp2=[false;temp1];

q(temp2)=0;

finalflag=zeros(size(newdates)); %1为买进信号，-1为卖出信号

finalflag(qix(q==1))=1;

finalflag(qix(q==-1))=-1;

buytime=qix(q==1);

soldtime=qix(q==-1);

buytimes=newdates(buytime); %买入时间点

soldtimes=newdates(soldtime); %卖出时间点

finalprice=bb(:,2);

buypoint=finalprice(finalflag==1); %买入价位

soldpoint=finalprice(finalflag==-1); %卖出价位

meanbuy=mean(buypoint);

meansold=mean(soldpoint);

legend('股票价格','20天平均线','葛兰碧法则1','葛兰碧法则3','葛兰碧法则4','葛兰碧法则5','葛兰碧法则7','葛兰碧法则8');

disp(['平均买入价格为',num2str(meanbuy),'元']);

disp(['平均卖出价格为',num2str(meansold),'元']);

D. Escape before bubbles crack

1. Introduction:

The bubble in the stock market cracked in July 2015, people were entrapped in the repentance for not escaping the market early. After all, we all took the bitter lesson in 2007 but people seemed to forget it when facing a boom. So I was curious. If people knew the bubble was going to collapse, how could they find a right junction to leave the stock market. So here I designed a BP neural network algorithm to learn the experience in 2007 and trid to apply those knowledge to 2015.

2. Codes:

a) Script 1:

%股票数据下载

close all

clear

load StockCodes.mat

tic

% n=size(Codes,1);

n=1032;

% a=ones(n,1)

% mat2cell(Codes,a)

flags=ones(n,1);%指示抓取数据是否成功

s=struct();

for i=1:n

try

c=yahoo;

stock=fetch(c,Codes(i,:),{'open','high','low','close','volume'},'2007/7/1','2007/12/1');

%雅虎抓取的数据是逆序的，就是最新的数据在最上面

dates=flipud(stock(:,1));

prices=flipud(stock(:,2:end));

s(i).prices=prices;

s(i).dates=dates;

disp(['抓取第 ' num2str(i) ' 只股票成功'])

catch %抓取失败，进行标记

flags(i)=0;

disp(['抓取第 ' num2str(i) ' 只股票失败'])

end

% i

end

toc %抓到852只股票，用时2900s

save sss3.mat s flags

function [ effects] = test_of_escape_BP_stock_fitness( perct )

%UNTITLED 此处显示有关此函数的摘要

% 此处显示详细说明

load BP_stock_escape2.mat

s_num=100;

result_exit=zeros(s_num,1);

quality_exit=zeros(s_num,1);

for i=(index(1:s_num))'

[ma,ma_index]=max(s(i).total(:,2));

[mi,mi_index]=min(s(i).total(:,2));

exit_index=find(s(i).probs>perct,1);

exit=s(i).total(exit_index,2);

if exit_index<ma_index

result_exit(i)=1;

quality_exit(i)=(exit-mi)/(ma-mi);

end

effects=-sum(quality_exit);

end

b) Script 2:

clear

load s_stock2015.mat; %存有07年股票数据 {'open','high','low','close','volume'}

load s_market2015.mat; %存有市场数据a

N=1032;

corrs=nan(N,1);

m_close=a(:,2);

m_temp=[false;diff(m_close)==0];

a(m_temp,2)=nan;

temp_fts1=fints(a);

temp_fts2=fillts(temp_fts1,'linear');

temp_fts3=toweekly(temp_fts2);

m_price=fts2mat(temp_fts3);

n2=size(m_price,1);

% ns=struct();

orders=find(flags);

nn=length(orders);

flags1=nan(N,1);

for j=1:nn

i=orders(j);

% volume=s(i).prices(:,5);

s_temp=s(i).prices; %载入股票各项变量

s_length=size(s_temp,1);

volume=s_temp(:,5); %交易量

temp=volume==0; %寻找无交易日

if sum(temp)<0.5*n2*5&s_length>0.5*n2*5; %筛掉数据较少的日子

s_temp(temp,:)=nan; %无交易日数据改为nan

% s(i).prices=s_temp;

temp_fts1=fints(s(i).dates,s_temp(:,4));

temp_fts2=fillts(temp_fts1,'linear');

temp_fts3=toweekly(temp_fts2);

s_price=fts2mat(temp_fts3);

n1=size(s_price,1);

n=min(n1,n2);

temp1=s_price(end-n+1:end);

temp2=m_price(end-n+1:end);

matrix=nancov(temp1,temp2);

[sigma,co]=cov2corr(matrix);

corrs(i)=co(1,2);

% s_temp(temp,:)=[];

s(i).prices(temp,:)=[];

disp(['第 ' num2str(i) ' 只股票处理成功'])

flags1(i)=i;

else

disp(['第 ' num2str(i) ' 只股票数据太少'])

end

a(m_temp,:)=[];

[new_corrs,index]=sort(corrs,'descend');

temp=[new_corrs,index];

temp1=isnan(new_corrs);

temp(temp1,:)=[];

index=temp(:,2);

s_num=100;

Input=[];

Output=[];

for i=(index(1:s_num))'

n1=20;

n2=3;

s_size=size(s(i).prices,1);

s_open=s(i).prices(:,1);

s_high=s(i).prices(:,2);

s_low=s(i).prices(:,3);

s_close=s(i).prices(:,4);

s_volume=s(i).prices(:,5);

for j=1:s_size-n1+1

s(i).price_20(j)=mean(s_close(j:j+n1-1));

s(i).volume_20(j)=mean(s_volume(j:j+n1-1));

s(i).dates_20(j)=s(i).dates(j+n1-1);

end

for j=1:s_size-n2+1

s(i).price_3(j)=mean(s_close(j:j+n2-1));

s(i).volume_3(j)=mean(s_volume(j:j+n2-1));

s(i).dates_3(j)=s(i).dates(j+n2-1);

end

% n=max(n1,n2);

trend=[s(i).price_20(1:end)',s(i).volume_20(1:end)',...

s(i).price_3(18:end)',s(i).volume_3(18:end)'];

trend=diff(log(trend));

temp1=size(trend,1);

[ma,ma_index]=max(s_close);

[mi,mi_index]=min(s_close(1:ma_index));

basis=s(i).prices(mi_index,:);

temp2=size(s(i).prices,1);

prices=s(i).prices./repmat(basis,temp2,1);

d=ma-mi;

f_v=1/2*sin((s_close-mi)/d)+1/2;

total=[s(i).dates(end-temp1+1:end),prices(end-temp1+1:end,:),trend...

f_v(end-temp1+1:end,:)];

s(i).total=total;

Input=[Input;total(:,2:end-1)];

Output=[Output;total(:,end)];

end

load BPnet_stock.mat

for i=(index(1:s_num))'

s(i).probs=(sim(net,s(i).total(:,2:end-1)'))';

end

save BP_stock_escape2.mat s index

load BP_stock_escape2.mat s index

test_of_escape_BP_stock_fitness( best_p )

c) Script 3:

s_num=100;

Input=[];

Output=[];

for i=(index(1:s_num))'

n1=20;

n2=3;

s_size=size(s(i).prices,1);

s_open=s(i).prices(:,1);

s_high=s(i).prices(:,2);

s_low=s(i).prices(:,3);

s_close=s(i).prices(:,4);

s_volume=s(i).prices(:,5);

for j=1:s_size-n1+1

s(i).price_20(j)=mean(s_close(j:j+n1-1));

s(i).volume_20(j)=mean(s_volume(j:j+n1-1));

s(i).dates_20(j)=s(i).dates(j+n1-1);

end

for j=1:s_size-n2+1

s(i).price_3(j)=mean(s_close(j:j+n2-1));

s(i).volume_3(j)=mean(s_volume(j:j+n2-1));

s(i).dates_3(j)=s(i).dates(j+n2-1);

end

% n=max(n1,n2);

trend=[s(i).price_20(1:end)',s(i).volume_20(1:end)',...

s(i).price_3(18:end)',s(i).volume_3(18:end)'];

trend=diff(log(trend));

temp1=size(trend,1);

[ma,ma_index]=max(s_close);

[mi,mi_index]=min(s_close(1:ma_index));

basis=s(i).prices(mi_index,:);

temp2=size(s(i).prices,1);

prices=s(i).prices./repmat(basis,temp2,1);

d=ma-mi;

f_v=1/2*sin((s_close-mi)/d)+1/2;

total=[s(i).dates(end-temp1+1:end),prices(end-temp1+1:end,:),trend...

f_v(end-temp1+1:end,:)];

s(i).total=total;

Input=[Input;total(:,2:end-1)];

Output=[Output;total(:,end)];

end

load BPnet_stock.mat

for i=(index(1:s_num))'

s(i).probs=(sim(net,s(i).total(:,2:end-1)'))';

end

save BP_stock_escape2.mat s index

load BP_stock_escape2.mat s index

test_of_escape_BP_stock_fitness( best_p )

[ma,ma_index]=max(s_close); %高点

[mi,mi_index]=min(s_close(1:ma_index)); %低点

basis=s(i).prices(mi_index,:);

temp2=size(s(i).prices,1);

prices=s(i).prices./repmat(basis,temp2,1); %标准化

d=ma-mi;

f_v=1/2*sin((s_close-mi)/d)+1/2; %逃顶函数值

total=[s(i).dates(end-temp1+1:end),prices(end-temp1+1:end,:),trend...

f_v(end-temp1+1:end,:)];

s(i).total=total;

Input=[Input;total(:,2:end-1)];

Output=[Output;total(:,end)];

end

net=newff(Input',Output',[9,5,1],{'purelin','logsig','purelin'});

net.trainParam.lr=0.05;

net.trainParam.epochs=200;

net.trainParam.goal=0.5*10^(-3);

net.divideFcn='';

net=train(net,Input',Output');

for i=(index(1:s_num))'

s(i).probs=(sim(net,s(i).total(:,2:end-1)'))';

end

save BP_stock_escape.mat s index s_num net

function [ effects] = escape_BP_stock_fitness( perct )

%UNTITLED 此处显示有关此函数的摘要

% 此处显示详细说明

load BP_stock_escape.mat

s_num=100;

result_exit=zeros(s_num,1); %逃顶结果

quality_exit=zeros(s_num,1); %逃顶质量

for i=(index(1:s_num))'

[ma,ma_index]=max(s(i).total(:,2));

[mi,mi_index]=min(s(i).total(:,2));

exit_index=find(s(i).probs>perct,1);

exit=s(i).total(exit_index,2);

if exit_index<ma_index

result_exit(i)=1;

quality_exit(i)=(exit-mi)/(ma-mi);

end

effects=-sum(quality_exit);

% options=gaoptimset('Generations',100,'Stallgenlimit',10,...

% 'plotfcn',@gaplotbestf);

% [x,f,reason]=ga(@escape_BP_stock_fitness,1,options);

perct=0.7:0.01:0.9;

effects=arrayfun(@escape_BP_stock_fitness,perct);

plot(perct,-effects);

xlabel('预警系数');

ylabel('预警效果');

[mi,mi_index]=min(effects);

-mi %最佳效果

best_p=perct(mi_index) %最佳系数 0.83

VI. Image Processing

I found the instruction Image Processing was so tempting that I learned it with one of friends. We work together to finish the required assignments in the class. It was exquisite, sutble and raletively pregmatic and I had a crush on it.

A. Pick the profile of corpuscles

1. Introduction:

This was the first assignment. We were asked to figure out the profile of the blood corpuscle which paved the road for further diagnosis. But professor noticed us that the former software was not capable to draw the cells in the edge of picture and told us he would be pretty happy if we could challenge this problem. And I did contrive an algorithm to achieve such an effect and I denominated it The Eraser.

2. Codes:

a=imread('blood512.bmp');imshow(a);

thresh=graythresh(a);

I=im2bw(a,thresh);

imshow(I);

m=edge(I,'canny');

m=double(m);m(m==0)=3;m(m==1)=0;m(m==3)=1;m=logical(m);%黑白颠倒；matlab对应的1为白色、0为黑色，故此时0为边框

%matlab对应的1为白色、0为黑色,即值越大，亮度越高。

[rows,cols]=size(m);

N=48; %最大的橡皮擦大小，即舍弃上限

step=2;%探测移动的步长

for n=16:4:N

imax=floor((cols-n)/step)+1;

jmax=floor((rows-n)/step)+1;

erases=ones(n);

erases(2:n-1,2:n-1)=0;

for i=1:imax

for j=1:jmax

x=(j-1)*step+1;%首点横坐标

y=(i-1)*step+1;%首点纵坐标

temp1=m(x:x+n-1,y);

temp2=m(x:x+n-1,y+n-1);

temp3=m(x,y:y+n-1)';

temp4=m(x+n-1,y:y+n-1)';

rim=[temp1;temp2;temp3;temp4];%边框之总

flag=any(rim==0);%边框是否被交叉，是返回1，否返回0

if ~flag

m(x+1:x+n-2,y+1:y+n-2)=1;

end

rim=[];

end

imshow(m)

B. Image Fusion

1. Introduction:

Professor gave a some exemplifications of former students in Image Fusion and he encouraged us to do it in different ways and so did I. The traditional way was the Wavelet Analysis, but here I conceived an Eigenvector method, composed of some knowledge in the class Linear Algebra.

2. Codes:

a) Script 1:

function NewPic=handlePic(filename,k)

%filename='Earth.jpg';

%k=10;

A=imread(filename);

A=double(A);

R=A(:,:,1);G=A(:,:,2);B=A(:,:,3);

%NewPic=zeros(size(A));

[u1,s1,v1]=svds(R,k);

[u2,s2,v2]=svds(G,k);

[u3,s3,v3]=svds(B,k);

im1=uint8(u1*s1*v1');

im2=uint8(u2*s2*v2');

im3=uint8(u3*s3*v3');

NewPic(:,:,1)=im1;

NewPic(:,:,2)=im2;

NewPic(:,:,3)=im3;

b) Script 2:

function NewPic=handlePic2(filename,k)

%filename='Earth.jpg';

%k=4;

A=imread(filename);

A=double(A);

R=A(:,:,1);G=A(:,:,2);B=A(:,:,3);

[c1,s1]=wavedec2(R,k,'sym5');

[c2,s2]=wavedec2(G,k,'sym5');

[c3,s3]=wavedec2(B,k,'sym5');

for i=1:k

Rx(:,:,i)=wrcoef2('a',c1,s1,'sym5',i);

Gx(:,:,i)=wrcoef2('a',c2,s2,'sym5',i);

Bx(:,:,i)=wrcoef2('a',c3,s3,'sym5',i);

end

for i=1:k

NewPic(:,:,1,i)=Rx(:,:,i);

NewPic(:,:,2,i)=Gx(:,:,i);

NewPic(:,:,3,i)=Rx(:,:,i);

end

%图像合成小波变换

a=imread('Earth.jpg');

b=imread('Ship.jpg');

n=4;

ax=rgb2gray(a);

bx=rgb2gray(b);

[c,s]=wavedec2(ax,n,'sym5');

for i=1:n

figure

a(:,:,i)=wrcoef2('a',c,s,'sym5',i);

imshow(a(:,:,i));

end

C. Counting Number

1. Introduction:

Our Professor wanted to create some software that can count the number of present students in the clas, which should be the grievous news for those truants. And we did achieve this goal but the accuracy was not high, which is the flaw of our design. I believe, however, we can fix this problem in the near future.

2. Codes:

a) Script 1:

a=imread('output_face.bmp'); %读取文件

se=strel('diamond',1); % 设置结构元素

temp1=imerode(a,se); %腐蚀

temp2=face_ferase(temp1,15,3); %橡皮擦

se=strel('diamond',2);

temp3=imdilate(temp2,se); %膨胀

%imopen()

%imclose()

[B,L] = bwboundaries(temp3,'noholes'); %连通域提取

n=0;

orders=[];

Limit=70;

for k = 1:length(B)

boundary = B{k};

if max(boundary(:,1))-min(boundary(:,1))>Limit...

||max(boundary(:,2))-min(boundary(:,2))>Limit

n=n+1;

orders=[orders;k];

L(L==k)=0;

end

imshow(label2rgb(L, @jet, [.5 .5 .5]))

hold on

for k = 1:length(B)

if ~any(orders==k)

boundary = B{k};

plot(boundary(:,2), boundary(:,1), 'w', 'LineWidth', 2)

end

N=size(B,1);

N-n %连通域的个数

figure

example001=B{16};

plot(example001(:,1),example001(:,2);

b) Script 2:

function m_new=face_ferase(b,N,step)

% N 最大的橡皮擦大小，即舍弃上限

% step 橡皮擦每次移动的步长，步长越小，精度越高，但运算速度越慢

%blood=imread('blood001.png');

%w=logical(blood);

[roww,colw]=size(b);

m=b;

%m=double(m);m(m==0)=3;m(m==1)=0;m(m==3)=1;

m=logical(m);%黑白颠倒；matlab对应的1为白色、0为黑色，故此时0为边框

m_new=m;

for n=15:5:N %15是边框初值

mm=true(roww+2*n,colw+2*n);

mm(n:n+roww-1,n:n+colw-1)=m_new;

m=mm;

%橡皮擦可以是各种各样的形状，四边形、六边形、八边形等等，这里仅使用四边形，边数越多、越接近于圆，效果越高

%step=2;%探测移动的步长

[rowm,colm]=size(m);

imax=floor((colm-n)/step)+1;

jmax=floor((rowm-n)/step)+1;

for i=1:imax

for j=1:jmax

x=(j-1)*step+1;%首点横坐标

y=(i-1)*step+1;%首点纵坐标

temp1=m(x:x+n-1,y);

temp2=m(x:x+n-1,y+n-1);

temp3=m(x,y:y+n-1)';

temp4=m(x+n-1,y:y+n-1)';

rim=[temp1;temp2;temp3;temp4];%边框之总

flag=all(rim==0);%边框是否包围住噪点

if flag

m(x+1:x+n-2,y+1:y+n-2)=0;

end

m_new=m(n:n+roww-1,n:n+colw-1);

end

imshow(m_new);

转载于:https://www.cnblogs.com/jessewei/p/5035769.html

weixin_30594001

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
Codes collection for my formal projects in university

In fact, these codes mean much more than a collection as they witnessed my growth.There is a better version, a delectable PDF brochure, which you can download in http://yun.baidu.com/share/li...
复制链接

扫一扫