Tomasz Mostowski

10.marca 2009

www.wne.uw.edu.pl/tmostowski

1

Matlab, zajęcia 3.

Pętle c.d.

Przypomnijmy sobie jak działa pętla for

Możemy podać normalnie w Matlabie

t=cputime;
for i=1:20

v(i)=i;

end
e=cputime-t

UWAGA: Taka operacja jest bardzo czasochłonna i nieoptymalna w Matlabie. W każdym
kroku pętli zmienia się wielkość wektora v. Tutaj służy wyłącznie do zaprezentowania
możliwości pętli.
Jeśli znamy długość wektora już przed operacją lepiej jest zrobić

clear v
t=cputime;
v=zeros(1,20);
for i=1:20

v(i)=i;

end

e=cputime-t

Zajmuje to zdecydowanie mniej czasu. Między innymi do tego służą macierze specjalne
(zeros, eye, ones). Proszę wykonać takie ćwiczenie dla macierzy o długości 500.

Jeszcze jeden przykład – metoda eliminacji Gaussa dla macierzy 3 na 3

B=[1 2 3; 4 10 12; 13 1 5]

%eliminacja w pierwszej kolumnie

for

j=2:3;

B(j,:)=B(j,:)-B(1,:)*B(j,1)/B(1,1);

end

%eliminacja w drugiej kolumnie

B(3,:)=B(3,:)-B(2,:)*B(3,2)/B(2,2);

B

Pętla While

Pętla while jest wykonywana tak długo jak spełniony jest jakiś warunek. Jest ona przydatna
szczególnie wtedy, kiedy nie wiemy jak długo będziemy wykonywać jakąś pętle.

Popatrzmy na przykład obliczeniowy

Tomasz Mostowski

10.marca 2009

www.wne.uw.edu.pl/tmostowski

2

Mamy macierz

1 1 1 1

0 0 0

0

0 0

0 0

0

0 0 0

e

A

e

e

e

⎡

⎤

⎢

⎥

⎢

⎥

⎢

⎥

=

⎢

⎥

⎢

⎥

⎢

⎥

⎣

⎦

Jeśli e=0, to macierz ma liniowo zależne kolumny i wtedy A’A jest macierzą osobliwą. Jeśli e
jest różne od zera macierz A’A jest odwracalna, gdyż A ma liniowo niezależne kolumny.
Przy obliczeniach numerycznych komputer dla małych e będzie zaokrąglał obliczenia tak
mocno, że A’A okaże się macierzą osobliwą.
Pytanie: dla jak małego e da się jeszcze obliczyć inv(A’A)?
Odp. Pętelka while

eps=1;

while

eps>1e-20

A=[ones(1,4);eps*eye(4)];

%stworzenie macierzy A poprzez połączenie

dwóch

eps

%pokaż gdzie jesteśmy?

d=det(A'*A)

%Wypisz wyznacznik. czy jest różny od zera?

eps=eps/10;

%zmniejsz eps

end

Można przejrzeć wyniki iż zobaczymy, że dla eps=1.0000e-008 wyznacznik już jest równy
zero.

Można jeszcze sprytniej.

eps=1;

while

eps>1e-20

A=[ones(1,4);eps*eye(4)];

%stworzenie macierzy A poprzez połączenie

dwóch

eps;

%pokaż gdzie jesteśmy?

d=det(A'*A);

%Wypisz wyznacznik. czy jest różny od zera?

if

d==0,

break

,

end

eps=eps/10;

%zmniejsz eps

end

eps

Proszę szczególnie zwrócić uwagę na linie (

if

d==0,

break

,

end

). Używamy tu opcji

break, która kończy wykonanie pętli while. Tak samo opcją break można kończyć pętle for.
W przypadku zagnieżdżonych pętli (pętli w pętli) break kończy tylko „swoją” pętle.

switch-case-otherwise

Jeżeli w zależności od wartości przyjmowanych przez jakąś zmienną program ma wykonać
różne czynności można użyć opcji switch-case-otherwise

Tomasz Mostowski

10.marca 2009

www.wne.uw.edu.pl/tmostowski

3

x=5;

switch

x

case

4

a=2*x;

case

5

a=3*x;

case

6

a=0;

otherwise

a=1;

end

a

Jeśli chcemy wykonać w kilku przypadkach jakąś czynność

x=6;

switch

x

case

{1,2,5}

a=3*x;

case

{6,4}

a=0;

otherwise

a=1;

end

a

Wczytywanie danych

W pracy ekonomisty zazwyczaj korzystamy z zewnętrznych zbiorów danych.
Popatrzmy jak można wczytać dane z plików excela.

num1 = xlsread('eur.xls')

Albo całą ścieżkę danych

num = xlsread('c:\dane\eur.xls');

Matlab wczytuje dane ignorując pierwszy wiersz.

Jeżeli chcemy dotrzeć do jakichś konkretnych informacji możemy po prostu wyciąć
odpowiednią kolumnę

close=num(:,4);

Możemy wczytać też tylko jakiś zakres danych

num = xlsread('eur.xls', 'C2:C100')

Notacja jest taka jak w Excelu.

Tomasz Mostowski

10.marca 2009

www.wne.uw.edu.pl/tmostowski

4

W ten łatwy sposób możemy wczytać dane liczbowe. Matlab ignoruje jednak dane tekstowe i
nie interpretuje polskiej daty z Excela. Bardzo często potrzebujemy jednak datę. Jak ją
wczytać? Należy otworzyć Excela i zmienić format daty na liczbowy. Następnie wczytać
dane do Matlaba. Np.

data = xlsread('eur.xls', 'B:B');

Kolejnym problemem jest inna interpretacja daty w Excelu i w Matlabie. Excel numeruje datę
od 1.stycznia 1900, natomiast Matlab od 1.stycznia 0000.

datestr(data(1))

Data jest wyświetlona jako 1.stycznia 0099.
Aby data była poprawnie wyświetlana należy „dodać” datę.

data=data+datenum('30DEC1899');

Sprawdźmy

datestr(data(1))

Teraz jest ok.


Wczytywanie danych w formacie .csv

Dane np. z portalu bossa.pl są zazwyczaj w formacie tekstowym oddzielonym przecinkiem,
czyli formatem CSV. Matlab ma funkcję csvread, ale wczytuje ona tylko dane numeryczne.

M = csvread('wig20.txt'); %błędne

Musimy podać numer wiersza i numer kolumny od której Matlab ma zacząć (W CSV
numeruje od zera).

M = csvread('wig20.txt',1,1);

Jak zwykle mamy problem z datą.

data=M(:,1);

datestr(data(1))

Nie wygląda to dobrze.
Żeby wczytać datę trzeba się trochę namęczyć np.

Najpierw zmieńmy liczbę na datę

Tomasz Mostowski

10.marca 2009

www.wne.uw.edu.pl/tmostowski

5

datastr=num2str(data);

rok=datastr(:,1:4);

mies=datastr(:,5:6);

dzien=datastr(:,7:8);

Żeby zmienić teraz w format daty można np. tak

sfdata=[mies(1,1:2),'/',dzien(1,1:2),'/',rok(1,1:4)]

datestr(sfdata)

Teraz jest odpowiednio. Jednym z zadań dzisiaj dla Państwa jest napisanie funkcji, która
wczytuje odpowiednio dane z zadanego formatu.

Można też inaczej

dataok=datenum(num2str(data),'yyyymmdd');
%zmiana danych na format tekstowy, zmiana na numeryczny używając zadanego formatu.
datestr(dataok(1))

Podstawy statystyki w Matlabie

W Matlabie mamy wbudowaną większość podstawowych statystyk
Popatrzmy na kilka podstawowych dotyczących np. zmiennej close

close=M(2:end,5);

/*opóźniona zmienna close*/
dclose=M(1:end-1,5);

/*stopa zwrotu, jako różnica logarytmów*/
zwrot=log(close)-log(dclose);

%średnia
mean(zwrot)

%wariancja
var(zwrot)

Dla zmiennej volume

volume=M(:,6);

Tomasz Mostowski

10.marca 2009

www.wne.uw.edu.pl/tmostowski

6

Odchylenie standardowe

std(volume)

%skośność
skewness(volume)
%kurtoza
kurtosis(volume)

3

1

3

4

1

4

(

) /

ˆ

(

) /

ˆ

N

i

i

N

i

i

x

x

N

skew

x

x

N

kurt

σ

σ

=

=

−

=

−

=

∑

∑

skewness(zwrot)
%kurtoza
kurtosis(zwrot)

O ile wiem wydziałowa wersja Matlaba nie ma Toolboxa zawierającego komendy takie jak
quantile czyli kwantyl.
Jak wiemy kwantyl rzędu p zmiennej losowej X to taka liczba k, że

P(X<k)=p

Np. median to kwantyl rzędu 0.5, czyli taka liczba med., że

P(X<med)=0.5

Nie znając rozkładu możemy estymować kwantyle używając kwantyli empirycznych.
Policzymy kwantyle dla dziennych stóp zwrotu. Dzienne stopy zwrotu obliczymy jako
różnicę logarytmów pomiędzy dwoma dniami. Na początek kilka prostych obliczeń.

[n m]=size(close)
%stworzenie zmiennej
dlclose=zeros(n,1);
%uzupełnienie jej liczbami
for i=2:n

dlclose(i)=log(close(i))-log(close(i-1));

end

Żeby obliczyć kwantyle empiryczne najłatwiej jest najpierw posortować dane

sr=sort(dlclose);

q25=sr(floor((n-1)/4))

Tomasz Mostowski

10.marca 2009

www.wne.uw.edu.pl/tmostowski

7

%Wartość znajdującą się w ¼ obserwacji

%3 kwartyl
q75=sr(floor(3*(n-1)/4))

Widzimy, że odległość między kwartylowa wynosi [-0.0099; 0.0105].

Tomasz Mostowski

10.marca 2009

www.wne.uw.edu.pl/tmostowski

8

Zadania

Zadanie 1.

Na podstawie kodu poznanego na zajęciach proszę napisać funkcję robiącą eliminację

Gaussa dla dowolnej macierzy kwadratowej (o dowolnych wymiarach).

Zadanie 2.
Napisać jeszcze jeden kod badający dla jak małego e macierz A’A nie jest
odwracalna. Użyj opcji while, ale uzależnij jej wykonywanie od wartości inv(A’A), a nie od
e.

1 1 1 1

0 0 0

0

0 0

0 0

0

0 0 0

e

A

e

e

e

⎡

⎤

⎢

⎥

⎢

⎥

⎢

⎥

=

⎢

⎥

⎢

⎥

⎢

⎥

⎣

⎦

Zadanie 3.
Wczytać dane z pliku bos.xls do Matlaba i stworzyć wektor danych kursu zamknięcia.

Zadanie 4.
Napisać funkcję wczytującą do Matlaba dane z pliku csv.

Zadanie 4a.
Wczytać dane z pliku TUP.txt

Zadanie 5.
Jeśli nie ma na wydziałowej licencji nie ma statystyk obliczających skośność i kurtozę
to stwórz funkcję obliczające je.

Ogólny raczej nie działający kod
h1=sum((ones(n,1)* mean(zmienna)-zmienna)^4)/n/std(zmienna)^4;

Zadanie 6.
Napisz

funkcję obliczającą kwantyl dowolnego rzędu. Niech argumentami funkcji

będzie zmienna i rząd kwantyla.

Zadanie 7.

Z pliku TUP.txt obliczyć dzienną stopę zwrotu (jako różnicę logaytmów kursu

zamknięcia) i obliczyć jej średnią, wariancję, kurtozę, skośność, medianę, kwantyle rzędu
0.05, 0.1, 0.15, 0.25, 0.75, 0.85, 0.9, 0.95