search

Charlotte324 / matlab-programming-and-numerical-method

1 stars 0 forks source link

matlab intro #3

Open Charlotte324 opened 1 year ago

Charlotte324 commented 1 year ago

black.m B-S 期权定价最基础的公式。

% Implements the Black pricing formula for a call option
% F: forward price for maturity T
% K: strike
% T: time to maturity
% r: risk free rate, continuosly compounded
% vol: volatility
function price = black (F, K, T, r, vol)
   B = exp(-r .* T);
   stdev = vol .* sqrt(T);
   d1 = log(F./K) ./ stdev + 0.5 .* stdev;
   d2 = d1 - stdev;
   nd1 = normcdf(d1);  % the function normcdf is available in Matlab
   nd2 = normcdf(d2);
   price = B .* ( F .* nd1 - K .* nd2 );
end
Charlotte324 commented 1 year ago

MATLAB

参考资料：

[mathworks get start with matlab](https://ww2.mathworks.cn/help/matlab/getting-started-with-matlab.html)

[Fabio bitbucket](https://bitbucket.org/fabiocannizzo/fe5116/src/master/)

intro

快捷

clear %清楚所有工作区储存的变量
pwd %告诉我们工作路径（当前文件夹的位置）
ls %告诉我当前文件夹里有哪些文件

edit xxxxx.m %打开或者创建某个当前文件夹下的文件-

基础运算

 +,-,/,*,^ %Basic Operators and precedences
.*,./ %最好写的时候都加个.，这样即使是矩阵也能没有bug地运行。
 ==,>,<,~= %logical operator 
 0 and 1 %boolean types
 &, &&, |, || %boolean operators
 help && %help

%%对输出结果进行简单的控制（在命令行用完了之后以后的都采取这个模式）
format long %小数位数多一点

%%回忆和调用之前输入过的命令
%按up（打游戏的上下左右小键盘的上），在输命令到一半的时候可以智能提醒

max();sqrt();log();exp();pi;normcdf()%需要下toolbox;
a == b
%当ab是矩阵的时候，只有当a和b矩阵行列数量一样的时候才能比（不然会输出矩阵维度必须一致error)
%当a,b的矩阵维度一致的时候，会生成与a,b同样大小的Logical值矩阵，里面是分别比较了a,b对应位置的每一个数是否相等。数只要相等，不管是不是0都会输出1的logical值。

&|~

我理解的是，返回与A的boolean相反的boolean。比如A是非零数，那么~A为0.

~(7>8)
%%%
1

要素

类型：Type-system (double, single, boolean, string, matrix)

keyword: type names; end; for;while 这些不能被用来作variable的名字

一般默认的数字都是double；

who %我有那些变量
whos %每个变量的信息

x=3 %bytes=8
x=single(3) %这样生成的x的type就不是默认的double而是single了。bytes=4
x=int64(3) %type=int64 bytes=8
x=int32(3) %type=int32 bytes=4
x=uint32(3) %type=uint32 bytes=4

1e-3 %1*10^(-3)
ok= sth==1
~ok =sth!=1 %由上面的默认。

符号

; %加在代码后面，可以防止运行的时候生成过多的结果（echo）

矩阵

建立矩阵，combine矩阵

[1,2;3,4]
[1 2;3 4]
[1 2
3 4]

ones(行,列) %全1矩阵
zeros(行,列) %全0矩阵
eye(维度) %单位矩阵
rand(1,2) %一行两列的random num

%%combine
z=[x x;y y] %只要串联的数组的维度一致就能combine
z=[zeros(4,4) eye(4);ones(1,8)]
%%%%%%%%%%%%%%%%%%%%%%%
ans =

     0     0     0     0     1     0     0     0
     0     0     0     0     0     1     0     0
     0     0     0     0     0     0     1     0
     0     0     0     0     0     0     0     1
     1     1     1     1     1     1     1     1
x='hello'
y='world'
m=[x ' ' y] %'hello world'

k='the value of sth is'
c=65 %对应的字母是A 66对应的字母是B
n=[k ' ' c+1] %'the value of sth is B'
o=[k ' ' num2str(c)]  %'the value of sth is 65'
q=[k ' ' 98] %'the value of sth is b' 因为98对应的是b
o=[k ' ' num2str(98)]

矩阵性质

size(x) %输出矩阵的行，列
size(x,1) %输出矩阵的行
size(x,2) %输出矩阵的列

矩阵变换

x
x' %转置
x+3 %x里每一个元素都加3
x*5 %x里每一个元素都乘5
x.^3 %x里每一个元素都变成自己的三次方

max(矩阵)

x是一个矩阵，max(x)可以输出每个column（列）的最大的那个值，生成一个一行的矩阵y。如果接着对y求max(y)，会对这行取一个最大的值。

max(a,b)，a和b是数的话会输出最大的那个数。

max(a,x)会输出矩阵z，z的每一个数都取max（a,原来x的对应位置的数）。matlab认为矩阵比数大。

x=[1 2;3 4; 5 6]
max(x,3)
%%%%%%%%%%%
ans =
     3     3
     3     4
     5     6

矩阵乘法

x*x %真正的乘法
x.*x %每一个元素乘以对应位置的元素 

inv(x) %矩阵的逆
x^(-1) %矩阵的逆

提取元素(element)

注意行和列都是从1开始数而不是从0

x=[1 2 0 0; 
   3 4 0 0; 
   1 1 1 1]
x(2,3) %矩阵（行，列）；matrix(row,column)
y=x(:,2) %拿出x的第二列
y=x(2:3,2:3) %拿出交点的四个值（2行3行和2列3列）
y=x(:,[2 2 2]) %拿出3次第2列，组成一个含三个列向量的行向量矩阵。
y=x(:,[1 3]) %拿出第1列和第3列

x(1:2,2:end) %拿出第2列到最后一列（第四列）的前两行
x(1:2, end:-1:2) %拿出第2列到最后一列（第四列）的前两行，但是换了个顺序
%%%%%%%%
     0     0     2
     0     0     4

赋值

提取元素=某值

x(1:2,[2 3 4])=1000 %把提取出来的6个元素的位置都换成1000
x(1:2,[1 3])=[10,20;30,40] %把提取出来的4个元素的位置换成右侧的矩阵对应位置的数字

等差数列矩阵

start num:差:end num会生成一行等差数列(row vector)（省略差的话默认为1）

end num不存在与等差数列的差的规律中时，就不会出现。

如果想要生成一列的矩阵(column vector)，直接加一个transpose（'）就好了。

x=1:3:11 %行向量的等差数列；超过11的不显示了
%%%%%%%%
     1     4     7    10
h=10:-1:1 %差也可以是负数
y=(1:10)' %列向量 ‘是转置符号

h=1:10
h(10:-1:1) 
h(end:-1:1)
%%%%%%%%
    10     9     8     7     6     5     4     3     2     1

解方程

AX=b 当成二元一次方程，A为系数矩阵，b是列向量，x是所求的俩变量的列向量。

A=[1 1;1 -1]; 
b=[2 0]';
x=A\b %[1 1]

显示

可以直接使用矩阵[]显示，也可以用fprintf和sprintf.

fprintf:[将数据写入到一个文本文件](。

sprintf:[将数据格式化为字符串或字符向量](https://ww2.mathworks.cn/help/matlab/ref/sprintf.html)。

fprintf(fileID,formatSpec,A1,...,An)

str = sprintf(formatSpec,A1,...,An)
[str,errmsg] = sprintf(formatSpec,A1,...,An)
str = sprintf(literalText)

例子：

%某 是说这个东西后面将用另一个东西来指代，这里就把它当成一个输入格式的地方就好了。

sprintf('%0.5g',1.61845624 ) %'1.6185' 从左往右留5个有效数字
sprintf('%0.5f',1.61845624 ) %'1.61846' 保留5位小数，4舍5入
sprintf('%70.5f',1/eps) %70位整数，5位小数。整数拿空格凑，小数拿0凑
%'                                                4503599627370496.00000'
sprintf('%s','hello') %'hello' 变成str
sprintf('The array is %dx%d',2,3) %'The array is 2x3'

[conversion character](https://ww2.mathworks.cn/help/matlab/matlab_prog/formatting-strings.html)

fprintf('the num is %d',63) 
%%%%%%%%
the num is 63>> 

fprintf('the num is %d\n',63) %\n可以加上换行
%%%%%%%%%%%%%
the num is 63
>> 

fprintf('the result of %d +%d is %d\n',2,3,2+3) 
%%%%%%%%%% 2+3会当成一个数自动给计算掉（因为最后一个的格式是d，数字而不是文本）
the result of 2 +3 is 5

mag=fprintf('the result of %d +%d is %d\n',2,3,2+3) 
%mag=24 字符数（空格和换行也算）

结构体 struct

里面的variables可以是不同种类

'a' 这种是label；

后面跟的那个是element。

s=struct('name','X1','value',10) %赋值之间用逗号隔开
s.value %10
s.name  %'X1'

s=struct('a',1,'b','sss',c=[1 2]) %我的运行不了，老师的这样可以
s=struct('a',1,'b','sss','c',[1 2]) %我的只有这样可以运行
%%%%%%%%%%%%
    a: 1
    b: 'sss'
    c: [1 2]

胞 cell

index+element

建立

a={1,[1,2,3,4]}
%%%%%%%
  1×2 cell 数组   %一行两列的意思
    {[1]}    {1×4 double}

调取

a{1,1} %1  拿出上面的1×2 cell 数组里面的第一行，第一列的东西
a{1,2} %[1,2,3,4]

函数

匿名函数

函数名字=@(x) 关于x的函数

positivepart=@(x) 3*x+1 %会在工作区生成一个function type的函数。
positivepart(2) %7
f=@(F)black(F,1,2,5/100,30/100) %black是同文件夹下的另一个定义的函数。这样可以嵌套着定义。

f=@(x,y) x.*x+y.*y;
f(2,3) %13

文件夹内文件写函数

function_name.m

%%%%%%%%%%%%%%在叫mysum的文件里
function y=mysum(a,b) 
    y=a+b+c; %最后不一定要输出y
end

%%%%%%%%%%%%%%在terminal里
mysum(1,2) %3

%%%%%%%%%%%%%%在另一个同文件夹下的xxx.m文件里
a=mysum(1,2) %a=3

一个有点综合的小例子：(画图)

%首先确认文件夹里有老师的black.m
function price = black (F, K, T, r, vol)
   B = exp(-r .* T);
   stdev = vol .* sqrt(T);
   d1 = log(F./K) ./ stdev + 0.5 .* stdev;
   d2 = d1 - stdev;
   nd1 = normcdf(d1);  % the function normcdf is available in Matlab
   nd2 = normcdf(d2);
   price = B .* ( F .* nd1 - K .* nd2 );
end

%%%%%%%%%%%%%%在当前.m里
f=@(F)black(F,1,2,5/100,30/100) 
F=[0.7:0.05:1.2];plot(F,f(F),F,max(0,F-1))

递归

recursive 定义函数时，也用到了函数本身

用阶乘举例：

function y=fact(n)
   assert(n>=0, 'n must be positive, but you entered: %d\n', n );
   assert(n == uint64(n), 'n must be integer, but you entered: %d\n', n );
   % 确认n可用了之后，我们开始递归
   y = fact2(n);
end

function y=fact2(n) %这个.m文件里有两个函数，但用到的是和.m同名的fact()。fact2只是为了方便定义fact()，在主文件里面没有办法使用fact2()
   if (n <= 1) %0和1的阶乘都是1
      y = 1;
   else
      y = n * fact2(n - 1);
   end
end

%为什么要使用两个函数：因为在第一个函数里有两个assert，如果把y=fact2(n)换成下面的if等的话，下面的if之类的因为用到了recursive所以变成了一个大循环。算fact(n-1)的时候，又会穿过前面的两个assert。这样写是比较节省算力的。

正常的定义方法（iterative），不用递归：

function f=fact(n)
    f=1
    for i=2:n %当n<2的时候，这个是empty box
        f=f*n
    end
end

画图

[画图](https://ww2.mathworks.cn/help/matlab/learn_matlab/plots.html)

plot(x,y) %x,y是一样长的行向量

x=0:0.1:2*pi;
plot(x,sin(x),x,cos(x)) %plot(横坐标，纵坐标，横坐标，纵坐标) 默认把第二个图和第一个图画在一个坐标系里

plot(x,sin(x),'.',x,cos(x),'g') %每个图后面加上这个图对应的格式要求。比如第一个图就是dot，点状的；第二个图是g，green绿色的

loglog(x,sin(x)) %两个都log化

语句结构

control flows

if else end

if ...  %不用冒号，用逗号或者直接不用任何符号都可以；condition用不用括号都可以
    ...
else   %也可以不用else
    ...
end

循环

loops(for while)

for i=1:10  
    do sth, %不用冒号，用逗号或者直接不用任何符号都可以。但写在一行的时候需要逗号
end

while (condtion)    %不用冒号，用逗号或者直接不用任何符号都可以；condition用不用括号都可以
    do sth
end

例子

for i=1:2:10,i,end
i=0;while (i<3),i,i=i+1,end %0，1，1，2，2，3 
i=0;while (i<3),i,i=i+1;end %0，1，2 %因为;的存在让结果不用输出出来

break continue

break：离开循环；

continue：终止现在的循环，去另一个循环。

for i=1:10,
    if(mod(i,2)==0)
        continue %回到for i=1:10 继续大循环
    end
    i
end
%%%%%%%%%%
1 3 5 7 9

for i=1:10,
    if (i>2)
        break %离开大循环，直接进入第二个end
    end
    i
end %第二个end

i=0;
while(1) %1作为一个非零数字永远是true的
    i,
    i=i+1;
    if(i>3)
        break %这个break导致了循环的结束
    end
end

debug

assert

是一个调试程序时经常使用的宏，在程序运行时它计算括号内的表达式，如果表达式为FALSE (0), 程序将报告错误，并终止执行。 如果表达式不为0，则继续执行后面的语句。 这个宏通常原来判断程序中是否出现了明显非法的数据，如果出现了终止程序以免导致严重后果，同时也便于查找错误。

使用递归的小例子：

function y=fact(n)
   if(n<0)
    sprintf('n must be positive, but you entered: %d\n', n );
    error('');
   end
   %assert(n>=0, 'n must be positive, but you entered: %d\n', n ); 
   %%这个assert与上面的if结构类似作用。
   assert(n == uint64(n), 'n must be integer, but you entered: %d\n', n );
   y = fact2(n); 
end

function y=fact2(n) 
   if (n <= 1) 
      y = 1;
   else
      y = n * fact2(n - 1);
   end
end

debug

点击代码行序号，可以设置断点（红色的），这样代码只会进行到这一步之前。

然后在terminal执行操作

fact(3)
>> fact(3)
4      y = fact2(n);
K>> n  %输入我们想查看的信息
n =
     3  %%%%说明这一步的时候n=3
%也可以直接在工作区看value的变动。

小绿箭头告诉我们目前运行到的位置，按F11或者按step in（步入）可以再往前前进一步。然后可以根据前进的过程来判断到底是哪里出了问题。

在函数调用堆栈里面，可以看到目前行进到了哪一块。随着用到的function不同，这里的内容也会发生改变（比如在没有运行到fact2()的时候，这个地方的fact2就会消失。

算力

tic,语句;toc

function dt=measureTime(fun, nRepeat, dryRun, args)
    % dry run: one-off initialization costs (e.g. compiling the m file)
    if (dryRun), fun(args); end

    tStart=tic;
    % total cost of the loop is nRepeat*costOf(fun)
    % if nRepeat is large enough, this will dominate on the cost of time()
    % and will be material with respect to the accuracy ot time()
    for i=1:nRepeat,
        fun(args);
    end
    dt=toc(tStart)/nRepeat;

其他小函数

sort(矩阵) %按顺序排