Chapter C. An Introduction to Matlab


Copyright 2016 Jon Danielsson. Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at. http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.

Listing C.1: Variables
Last edited: 2011

x = 10;
y = 5;
z = x+y;
who
Your variables are:
ans  x    y    z
		

Listing C.2: Arrays
Last edited: 2011

a = 1:5
a =
  1     2     3     4     5
b = [1,2,3,4;5,6,7,8]
b =
     1     2     3     4
     5     6     7     8
size(b)
ans =
     2     4
		

Listing C.3: Sequences
Last edited: 2011

x = 0:2:10
x =
     0     2     4     6     8    10
x = linspace(0, 10, 6)
x =
     0     2     4     6     8    10
x = logspace(0,2,4)
x =
    1.0000    4.6416   21.5443  100.0000
		

Listing C.4: Load files
Last edited: 2011

sp=load('sp.dat'); 
		

Listing C.5: Download S&P-500 data in Matlab
Last edited: 2011

price = hist_stock_data('01012000','31122000','^gspc');  
y=diff(log(price.Close(end:-1:1)))  
plot(y)                
		

Listing C.6: Basic data analysis
Last edited: 2011

length(y)        
sum(y)          
prod(y)          
range(y)        
mean(y)          
median(y)        
var(y)          
std(y)          
sqrt(var(y))    
corrcoef(y)      
skewness(y)      
kurtosis(y)      
quantile(y,0.01)    
min(y)          
max(y)          
sort(y)          
abs(y)          
diff(y)          
		

Listing C.7: Matrix calculations
Last edited: 2011

A = [2,4;7,3]
B = [3,2;1,4]
A .* B
ans =
     6     8
     7    12
A * B
ans =
    10    20
    24    26
C = [1,2]
A * C  
??? Error using ==> mtimes
Inner matrix dimensions must agree.
A * C'    
ans =
    10
    13
diag(C)    
ans =
       1     0
       0     2
eye(3)    
ans =
     1     0     0
     0     1     0
     0     0     1
		

Listing C.8: Useful matrix functions
Last edited: 2011

M = [1 2;3 6]    
rank(M)          
inv(M)          
det(M)          
eig(M)          
[V,D] = eig(M)  
sqrtm(M)        
trace(M)        
null(M)          
poly(M)          
chol(M)          
		

Listing C.9: Distribution functions
Last edited: 2011

q = 0                  
p = 0.5                
mu = 1                
sigma = 2             
normpdf(q,mu,sigma)   
normcdf(q,mu,sigma)   
norminv(p,mu,sigma)   
normrnd(mu,sigma,10,2) 
normfit(x) 
		

Listing C.10: Common distributions
Last edited: 2011

trnd(3)                 
chi2rnd(3)              
unifrnd(1,2)            
lognrnd(1,2)        
exprnd(2)              
gevrnd(10,1,2)      
gprnd(10,1,2)    
binornd(3,0.3)            
poissrnd(0.4)
		

Listing C.11: Multivariate distributions
Last edited: 2011

mvnrnd(mu,sigma) 
mvtrnd(C,df)     
copularnd('Gaussian',rho,N)
copularnd('t',rho,NU,N)
copularnd('Clayton',alpha,N)
		

Listing C.12: Testing for normality
Last edited: 2011

alpha=0.01
[h,p,jbstat,critval] = jbtest(y,alpha)   
h     
p     
jbstat = 20.0648  
critval = 11.5875     
[h,p,ksstat,critval]=kstest(y,[],alpha)  
h = 1
p = 0
ksstat = 0.4814
critval = 0.1020
qqplot(y)   
qqplot(x,y) 
		

Listing C.13: Basic time series analysis
Last edited: 2011

[ac,acstd] = sacf(y,20,1,0);       
		

Listing C.14: A script M-file
Last edited: 2011

a = randn(10);
b = a * 5 + 2;
plot(b)
		

Listing C.15: A function M-file
Last edited: 2011

function  k = mykurtosis(x)
  m4 = mean((x - mean(x)).^4);
  k = m4/std(x)^4 - 3
end  
		

Listing C.16: Calling function M-file
Last edited: 2011

mykurtosis(y)
k = 1.3502
		

Listing C.17: Kurtosis function
Last edited: 2011

function  k = mykurtosis1(x,a)
  if nargin == 1
    a = 3;
    end
  m4 = mean((x - mean(x)).^4);
  k = m4/std(x)^4 - a
		

Listing C.18: A for loop
Last edited: 2011

x = randn(1,5);
z = NaN(1,5);  
for i=1:5
  z(i)=x(i)+0.5;
end
		

Listing C.19: An if-else loop
Last edited: 2011

a = 10;
if (rem(a,3))==0
  disp('a is a multiple of 3')
else
  disp('a is not a multiple of 3')
end  
		

Listing C.20: A while loop
Last edited: 2011

a = 1;
n = 1;
while a < 100
  a = a*n;
  n = n + 1;
end
		

Listing C.21: A switch-case loop
Last edited: 2011

x = 30;
units='kg';
switch units
  case{'kilograms','kg'}
      y = x*1000;
  case{'grams','g'}
      y = x;
  case{'ounce','oz'}
      y = x*28.35;
  case{'pounds','lb'}
      y = x*453.6;
  otherwise
      disp(['Unknown Units:' units])
      y = NaN;
end
		

Listing C.22: Loop avoidance
Last edited: 2011

a = 1:1000
b = 1000 - a
ssum = 0;
for i = 1:1000 
   ssum = ssum + a(i)*b(i);
end
ssum = a*b'
		

Listing C.23: Maximum likelihood estimation
Last edited: 2011

randn('seed',1);   
N = 100;            
x = randn(N,1)*2+3;  
theta0 = [-2, 5];    
global x;            
[theta, likelihood_value] = fminunc(@likfunc, theta0) 
theta = 2.8082    1.9366
likelihood_value = 116.0928
		

Listing C.24: Normal log likelihood
Last edited: 2011

function loglik = likfunc(theta)
  global x;
  T = length(x);
  mu = theta(1);
  sigma2 = theta(2)^2;
  loglik = 0.5*T*log(sigma2);
  loglik = loglik +  0.5*(sum((x-mu).^2 / sigma2));