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