Z Skrypty dla studentów Ekonofizyki UPGOW

Wersja z 16:15, 13 lut 2010

Analiza Szeregów Czasowych

Definicja szeregu czasowego

Możemy spotkać różne definicje szeregu czasowego.

Szereg czasowy to

• ciąg obserwacji pokazujący kształtowanie się badanego zjawiska w kolejnych okresach czasu (sekundach, dniach, latach, itp.).
• realizacja procesu stochastycznego, którego dziedziną jest czas; to ciąg informacji uporządkowanych w czasie, których pomiary wykonywane są z dokładnym krokiem czasowym.
• ciąg obserwacji x_t zapisywanych w ściśle określonym czasie.

Wśród składników szeregu czasowego możemy wyróżnić:

• trend (tendencję rozwojową),
• wahania sezonowe,
• wahania cykliczne (koniunkturalne),
• wahania przypadkowe.

W jakim celu badamy szeregi czasowe?

Analiza tego typu zagadnień ma generalnie dwa podstawowe cele:

1. odgadnięcie natury danego zjawiska losowego, tj. badanie własności szeregu i znalezienie modelu najlepiej opisującego zjawisko,
2. prognozowanie (predykcja), tj. przewidywanie kolejnych wartości szeregu czasowego na podstawie znalezionego modelu.

Przykłady szeregów czasowych

Przykład 1. Prąd płynący przez opornik.

Rysunek 1. 100 kolejnych punktów czasowych dla szeregu czasowego z przykładu 1 dla dwóch wartości oporu r.

Jeżeli do opornika charakteryzującego się oporem \(r\) przyłożymy zmienne napięcie

\( U(t) = a \cos (\omega t), \! \)

gdzie \(a\) to amplituda zmiennego napięcia przyłożonego do opornika, a okres zmienności to \(T = 2 \pi / \omega\). Wtedy natężenie prądu elektrycznego płynącego przez opornik można wyrazić wzorem

\( I(t) = \frac{a \cos (\omega t)}{r}. \! \)

Jest to oczywiście ciągła funkcja czasu, jednak, kiedy będziemy rejestrować wartości natężenia \(I(t)\) w kolejnych chwilach czasu (np. co \(0.1 T\), 1 milisekundę czy 1 godzinę), dostaniemy dyskretny szereg czasowy \(I_i\) indeksowany kolejnymi pomiarami \( i = 0, 1, 2, \dots \). Przykładowe szeregi czasowe opisane powyższym wzorem można znaleźć na rysunku 1.

Ćwiczenie W1.1

Wygeneruj w programie Matlab/Octave rysunek 1 (legenda jest opcjonalna).

1. Zbierz do tablic indeksy \(i\) oraz wartości natężenia prądu w punktach \(t_i = i \cdot ( 6 \pi / 100 ), i \in [0,100]\).
2. Wyplotuj do pliku (np: rysW11.png) wykres \(I_i = a \cos(\omega t_i + \phi) / r \).

close all
 
h = figure;
set (h,'papertype', 'a4')
set (h,'paperunits','centimeters');
set (h,'papersize',[8 6])
set (h,'paperposition', [0,0,[8 6]])
set (h,'defaultaxesposition', [0.15, 0.15, 0.75, 0.75])
set (0,'defaultaxesfontsize', 24)
 
x = -pi:7*pi/100:6*pi;
i = 0:1:100;
 
y = sin(x)/0.8;
plot (i,y ,"+^; r = 0.9 {/Symbol o}, {/Symbol f} = 0;")
 
hold on
 
y = sin(x + pi/3)/1.5;
plot(i,y,"-or; r = 1.5{/Symbol O}, {/Symbol f}= {/Symbol p}/3;");
 
xlabel('i');
ylabel('I_i');
title('I_i = a cos (\omega t_i + \phi) /r,  a = 1V, \omega = 1Hz'); 
 
grid on
 
print('example01m.eps','-deps','-FTimes-Roman:24');
print('example01m.svg','-dsvg');
print('example01m.png','-dpng');
%system('convert -density 100 example01m.eps example01m.png')

import matplotlib.pyplot as plt
import numpy as np
 
plt.figure(1, figsize=(8,6), dpi=600)
i = np.arange(100)
x = np.arange(-np.pi, 6.0*np.pi, 7.0*np.pi/100)
 
kolor  = ['b','r']
lines  = ['v','-p']
lw = 1
 
j=0
a, o, r, phi   = 1, 1, 0.8, 0
labels = r'$r = %.1f\Omega, \quad \phi = %.1f$' % (r,phi)
y = a * cos(o*x + phi) / r
plt.plot(i, y, lines[j], color=kolor[j], label=labels, linewidth=lw)
 
j=1
a, o, r, phi   = 1, 1, 1.5, np.pi/3.
labels = r'$r = %.1f\Omega, \quad \phi = \pi/3$' % (r)
y = a * cos(o*x + phi) / r
plt.plot(i, y, lines[j], color=kolor[j], label=labels, linewidth=lw)
 
plt.xlabel(r'$i$')
plt.ylabel(r'$I_i$')
plt.grid(True)
 
leg = plt.legend(shadow=True,
    fancybox=True,
    mode='expand',
    ncol=2,
    loc=(0.025,0.96),
    handletextpad=0.1, 
    title=r'$I_t = a \cos (\omega t + \phi) / r \qquad a = 1V, \quad\omega = 1 Hz$')
plt.savefig('example01.png')

Przykład 2. Proces dwustanowy (proces dychotomiczny, binarny, zerojedynkowy).

Rysunek 2. 100 kolejnych punktów czasowych dla rzutu monetą (przykład 2). Dane użyte do wygenerowania wykresu: 1 1 1 0 0 0 0 1 1 0 0 1 1 1 1 1 0 1 0 1 1 0 0 0 0 1 1 0 0 1 1 1 1 0 0 1 0 0 0 1 0 1 0 0 1 0 1 0 1 0 1 0 1 1 0 1 1 0 1 1 0 0 0 1 0 0 1 1 0 1 0 1 1 0 1 0 0 0 0 0 1 1 0 1 1 1 0 1 0 0 1 0 1 1 1 0 1 0 1 1.

Niech \(\{X_t, t = 1,2,3,\dots\}\) będzie uporządkowanym zbiorem niezależnych zmiennych losowych (sekwencją losową), dla których prawdopodobieństwo

\( P (X_t = 0) = P (X_t = 1) = 1/2. \)

(dowód istnienia potrzebnej przestrzeni probabilistycznej na razie sobie darujemy). Seria pomiarowa składać się będzie z losowo ułożonych w czasie zer i jedynek {0,0,0,1,0,1,1,1,1,0,...}. Przykładem jest rzut monetą.

Ćwiczenie W1.2

Każda osoba ma za zadanie wykonać N (w zależności od liczebności grupy, w sumie około 100 na wszystkich studentów) rzutów monetą. W arkuszu kalkulacyjnym na [docs.google.com] wpisujemy wartości:

• 0 jeżeli wyrzuciliśmy Orła
• 1 jeżeli wyrzuciliśmy Reszkę

każdy w oddzielnej kolumnie. Stwórz prosty wykres danych w arkuszu kalkulacyjnym Google. Następnie za pomocą programu Matlab/Octave stwórz rysunek przedstawiający tak utworzony szereg czasowy. Szereg ma uwzględniać pomiary wszystkich.

• wyeksportuj dane z arkusza Google do pliku CSV
• zaimportuj dane do tabeli w Matlab'ie
• Wyplotuj do pliku (np: rysW12.png) wykres \(X_t\).

close all
 
h = figure;
set (h,'papertype', 'a4')
set (h,'paperunits','centimeters');
set (h,'papersize',[8 6])
set (h,'paperposition', [0,0,[8 6]])
set (h,'defaultaxesposition', [0.15, 0.15, 0.75, 0.75])
set (0,'defaultaxesfontsize', 24)
 
N = 100;
i = 0:N;
 
%
% eksperyment studentów (symulcja rzutu monetą)
% odkomentowac jezeli nie mamy danych
%
%y = int32(rand(N+1,1));
%csvwrite('dataW12.csv',y);
%
 
y = csvread('dataW12.csv');
plot (i,y ,"-^; Orzeł czy Reszka?;")
axis([0,N,-0.5,1.5]);
xlabel('t');
ylabel('X_t');
grid on
 
print('example02m.eps','-deps','-FTimes-Roman:24')

# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
from random import randint
 
plt.figure(1, figsize=(8,6), dpi=600)
i = np.array([s for s in np.arange(100)])
y = np.array([randint(0,1) for s in np.arange(100)])
 
plt.plot(i, y, 'v', color='b', label=u"Orzeł czy Reszka?", linewidth=1)
plt.axis([0, len(y), -0.5, 1.5])
plt.xlabel(r'$t$')
plt.ylabel(r'$X_t$')
plt.grid(True)
 
leg = plt.legend(shadow=True
    ,fancybox=True
    ,ncol=1
    ,loc=(0.025,0.96)
    )
plt.savefig('example02.png')

Przykład 3. Populacja Polski.

Rysunek 3. Zmienność populacji Polski w latach 1960 - 2010 (przykład 3). Dane: http://epp.eurostat.ec.europa.eu

Przykład zmiany liczby ludności Polski.

PierwszyRok = 1960;
OstatniRok  = 2010;
%N = OstatniRok - PierwszyRok;
%i = 0:N;
% lub
i = PierwszyRok:OstatniRok;
y = csvread('PopulacjaPolski1960-2010.csv');
 
plot(i,y(:,5:5)/1000.,'--r; Populacja Polski w latach 1960-2010;',
                'marker','^', 'MarkerSize',14, 'markeredgecolor','black')
 
xlabel('t');
ylabel('X_t (w tysiącach)');
grid on;
legend('Location','SouthEast');

TBA

Ćwiczenie W1.3: Populacja Polski w latach 1960 - 2010.

Eurostat’s mission is to provide the European Union with a high-quality statistical information service.

1. Strona Eurostat-u może posłużyć Państwu jako doskonałe źródło ciekawych danych statystycznych. Ze strony Eurostat-u proszę pobrać interesujące nas dane, tj. wygenerować plik *.csv zawierający dane dotyczące stanu liczebnego Polski z w latach 1960 - 2010. Z pliku zawierającego dużo więcej danych proszę wyodrębnić te właściwe i wyplotować do pliku.

1. Proszę w internecie poszukać jak najdalej wstecz sięgających danych statystycznych odnośnie ludności Polski i powtórzyć procedurę z tymi danymi (da się znaleźć dane od roku około 1000).

Źródło danych

• http://epp.eurostat.ec.europa.eu/portal/page/portal/statistics/search_database
• http://pl.wikipedia.org/wiki/Ludność_Polski

Przykład 4. S&P 500

Rysunek 4. Indeks S&P 500 z lat 1987 - 2010.

S&P 500 jest indeksem w skład którego wchodzi 500 firm o największej kapitalizacji, notowanych na New York Stock Exchange i NASDAQ, są to głównie firmy amerykańskie. Indeks ten jest najbardziej znanym wskaźnikiem zarządzanym przez Standard & Poor's (oddział McGraw-Hill). S&P 500 wchodzi w skład szerszego indeksu - S&P 1500 oraz S&P Global 1200.

• Lista firm należących do S&P 500 z procentowym udziałem w indeksie

Ćwiczenie W1.4: S&P 500.

1. Utwórz rysunek 4.

y = csvread('SandP500.csv');
 
% wszystkie dane na raz
% plot(y(2:size(y,1),7:7),'--r; S&P 500 1987 - 2010 (close);')
 
bN = 2;
N  = 165 ;
i = bN:N;
plot(i,y(bN:N,7:7),'--b');
hold on
 
bN = N;
N  = 196 ;
i = bN:N;
plot(i,y(bN:N,7:7),'--r');
hold on
 
bN = N;
N  = 250;
i = bN:N;
plot(i,y(bN:N,7:7),'--b');
 
 
xlabel('t');
ylabel('X_t');
grid on;
title('S&P 500 1987 - 2010 (close)');

# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
 
plt.figure(1, figsize=(8,6), dpi=600)
 
fd = open('SandP500.csv', 'rb')
y = np.loadtxt(fd, delimiter = ',', usecols=[6], unpack=True, skiprows=1)
i = np.array([i for i in range(0,len(y))])
 
# poniższe punkty pochodzą z wstępnej, wzrokowej analizy wykresu
hobe = [0,164,195,249,len(y)-11,len(y)]
 
kolorki = ['b','r']
legenda = ['hossa','bessa']
 
for hb in range(0,len(hobe)-1):
  if hb > 1:
    legenda = [None,None]
  plt.plot(i[hobe[hb]:hobe[hb+1]],y[hobe[hb]:hobe[hb+1]], '-^', color=kolorki[hb%2], label=legenda[hb%2], linewidth=1)
 
plt.xlabel(r'$t$')
plt.ylabel(r'$X_t$')
plt.grid(True)
 
leg = plt.legend(shadow=True
    ,fancybox=True
    ,ncol=2
    ,loc=(0.025,0.98)
    ,title=r'S&P 500, 1987 - 2010 (close)'
    )
plt.savefig('example04py.png')

Źródło

http://finance.aol.com/quotes/sandp-500-index-rth/$inx/cmi/historical-prices?tf=all&gran=d

Klasyczne przykłady z książki Brockwell-a

Liczba ludności USA, lata 1790 - 1980

Rysunek 5. Liczba ludności USA w latach 1790 - 1980, w odstępach dziesięcioletnich.

Miesięczna liczba śmierci w wypadkach, USA, lata 1973 - 1978