Что такое R?

R - интерпретируемый язык программирования, который был разработан специально для статистического анализа и анализа данных. Кроме этого важно знать, что R распространяется в рамках проекта GNU (то есть с открытым исходным кодом). Зачем он нам нужен:

  1. R очень легок в освоении.
  2. Пожалуй, самая красивая визуализация данных.
  3. Легкая статистическая обработка данных.
  4. Множество готовых пакетов, которые уже разработаны для решения вашей задачи.
  5. Большое комьюнити.
  6. R - это полнофункциональный язык программирования, который используется в различных областях науки: математика, биология, лингвистика, история, экономика, etc. В принципе некорректно сравнивать R с пакетами вроде gretl.

Да, если кто-то решил, что R - очередная фикция, то вот список компаний, которые используют R.

Интерпретируемый язык

Интерпретируемый язык программирования - такой, где команды исполняются одна за одной, что делает их использование более удобными и приятными. Самый простой вариант использования R - калькулятор. Давайте проведем пару математических расчетов. Для этого нам понадобится знать всего один самый главный оператор – оператор присваивания <-, он записывает в переменную некоторое значение. Ниже представлен небольшой фрагмент кода, который Вы можете исполнить на своем компьютере. После фрагмента идет результат выполнения этих команд.

a <- 42
a
## [1] 42

Здесь мы первой командой устанавливаем значение переменной a, а второй командой выводим ее значение. R напечатал не просто 42, а [1] 42. Пусть сейчас это вас не смущает. [1] – это намек на то, что полученное в результате работы значение - скаляр, а не вектор. Теперь создадим еще пару переменных и что-нибудь подсчитаем.

b <- 7
a*b
## [1] 294
cos(a)*log(a+b) # Можно считать и более сложные вещи :)
## [1] -1.556671

И так далее. Кстати, в предыдущем блоке кода есть комментарий. Все, что идет после символа # не будет выполняться. Комментировать свой код - правильно и хорошо.

Визулизация данных

В R очень просто визуализировать все что угодно. Просто не обращайте внимания на команды, если не понимаете их:

plot(sin(seq(0,2*pi,0.1)), col=2, cex=2, main="Синус" )

boxplot(rnorm(100, 0, 1), rnorm(100,1,5), main="Две выборки: из N(0,1) и из N(1,5)")

hist(rchisq(1000, 10), breaks = 30, col=3, main="Хи-квадрат с 10 степенями свободы")

Это только базовые примеры визуализации. А вот один из моих любимых примеров: визуализация gps-треков вело поездок по Лондону. К сожалению, не моих. Если кто-то сможет построить хотя бы в половину настолько же красивый график в Excel/Matlab/gretl/Maple/…, то с меня дополнительные баллы.

Пакеты для всего

Но настоящая сила R кроется в пакетах. Пакеты - готовые наборы функций для решения различных задач. В настоящее время в официальном репозитории пакетов CRAN хранится 7178 пакетов, разработанных для решения задач всевозможных областей.

Единственный минус – очень часто так и хочется просто взять готовый пакет, чтобы не напрягаться самому. В нашем курсе все ключевые моменты мы будем делать сами :)

Базовые команды

Пока что мы умеем только присваивать переменным скалярные значения (числа). Но базовым типом переменной в R является вектор - упорядоченная последовательность чисел. Вектор можно задать разными способами. Основные - это оператор : который создает последовательность целых чисел и команда seq, которая создает более гибкие вектора.

a <- 1:10 # Создает вектор чисел от 1 до 10
a
##  [1]  1  2  3  4  5  6  7  8  9 10
b <- seq(3, 4, by=0.05) # Вектор от 3 до 4, разделенный на отрезки по 0.05
b
##  [1] 3.00 3.05 3.10 3.15 3.20 3.25 3.30 3.35 3.40 3.45 3.50 3.55 3.60 3.65
## [15] 3.70 3.75 3.80 3.85 3.90 3.95 4.00

Теперь кстати понятно стало обозначение [1] в начале каждой строки. В последней команде вывод занял больше одной строки, и 15 в скобках указывает, что 15-ый элемент вектора – 3.70.

Еще два способа получить вектор - это команда rep и универсальная команда c:

a <- rep(0, times=10)
a
##  [1] 0 0 0 0 0 0 0 0 0 0
vect <- c(3, 6, 9, -2) # Создадим вектор из четырех чисел ...
vect
## [1]  3  6  9 -2
c(vect, 5, 7) # ... и допишем к нему еще два.
## [1]  3  6  9 -2  5  7

Все операции R старается выполнить векторно, поэлементно. Тут все должно быть ясно из примеров:

a <- 1:10 # Вектор от 1 до 10
b <- 1:10 # и еще один такой же
a+b
##  [1]  2  4  6  8 10 12 14 16 18 20
a-b
##  [1] 0 0 0 0 0 0 0 0 0 0
a*b
##  [1]   1   4   9  16  25  36  49  64  81 100
a^b
##  [1]           1           4          27         256        3125
##  [6]       46656      823543    16777216   387420489 10000000000
a+5
##  [1]  6  7  8  9 10 11 12 13 14 15
log(a)
##  [1] 0.0000000 0.6931472 1.0986123 1.3862944 1.6094379 1.7917595 1.9459101
##  [8] 2.0794415 2.1972246 2.3025851

А теперь время узнать про главную команду в R – ?. Да, это вопросительный знак. Эта функция вызывает документацию по любой другой функции, где всегда можно что-то уточнить и проверить. Наберите в своей рабочей среде:

?plot

Чтобы сейчас не закапываться в изучение базовых команд мы будем изучать их по ходу дела, потому что иначе они не усвоятся. И пока плавно перейдем к реальным задачам.

Анализ гистограм

Дальше в курсе мы почти все время будем анализировать распределение доходностей финансовых активов. Поэтому поговорим немного о распределениях и о анализе этих самых распределений, то есть гистограмм.

Немного о распределениях

Так как R был создан статистиками для статистиков, то все необходимые распределения уже есть в R. Необходимые – это: равномерное, нормальное, t (Стьюдента), Хи-квадрат, F-распределение, гамма, и многие другие. Каждое из этих распределений имеет в R свое название:

  • Равномерное – unif
  • Нормальное – norm
  • Стьюдента – t
  • Хи-квадрат – chisq

Для каждого из этих распределений существует несколько функций: r*, p*, d*, q*. Где звездочка – это имя распределения. Функции r* – это генераторы случайных чисел. То есть runif генерирует равномерно распределенное случайное число. А rnorm - случайное нормальное.

runif(3) # Единственный параметр - количество чисел, которые нужно сгенерировать.
## [1] 0.71414881 0.07147234 0.43141084

Параметры распределения необходимо также передавать функциям, иначе будут использоваться стандартные. Чтобы узнать какие параметры ожидает каждое распределение, не стесняемся открывать справку:

?rnorm

Хотим 14 нормальных чисел со средним 7 и стандартным отклонением 0.5? Так и пишем:

rnorm(14, mean=7, sd=0.5)
##  [1] 7.359320 6.419975 7.683645 5.679030 7.089875 7.517579 6.520223
##  [8] 6.370759 7.483771 6.859768 7.451023 6.932730 7.715658 7.122656

Функции q* (qnorm, qunif, …) возвращают квантили распределения. Какой там квантиль t-распределения с большим количеством степеней свободы на уровне 0.95?

qt(0.95, Inf) # Тут допустимо особое значение - Inf, бесконечное число
## [1] 1.644854

Функции p* и d* возвращают функцию плотности распределения и функцию плотности соответственно. Подробнее читайте в справке.

Анализируем данные

В R есть много встроенных датасетов (наборов данных). Возьмем для нашего просто анализа один из них - вес куриц в зависимости от типа корма. Но эту зависимость мы пока изучать не будем, а просто проанализируем распределение весов куриц.

w <- chickwts$weight # Сохраняем веста в переменную w'
hist(w, breaks = 20) # Строим гистограмму с 20 колонками

А команда ecdf умеет строить функцию плотности:

plot(ecdf(w), do.points=FALSE, verticals = TRUE) # Вы вольны изменять параметры

Найдем базовые характеристики распределения:

mean(w)
## [1] 261.3099
median(w)
## [1] 258
var(w)
## [1] 6095.503
sd(w)
## [1] 78.0737
quantile(w, probs = c(0.05, 0.5, 0.95))
##    5%   50%   95% 
## 140.5 258.0 385.0

Медиана близка к среднему: какой вывод можно сделать?

Схожесть распределений

Пусть нам показалось, что распределение весов куриц - нормальное. Как это проверять? Первый вариант - визуально:

hist(w, breaks = 20, freq = FALSE)
points <- seq(min(w), max(w), length.out = 100)
lines(points, dnorm(points, mean = mean(w), sd = sd(w)), col=2)

Близко ли распределение к нормальному?

Второй и более солидный вариант – статистические тесты. Во-первых, тест Шапиро-Уилка. Гипотеза \(H_0: x \sim N(\mu, \sigma)\). Статистику можно глянуть в сети. В R для проведения теста само собой есть готовая функция.

shapiro.test(w)
## 
##  Shapiro-Wilk normality test
## 
## data:  w
## W = 0.97674, p-value = 0.2101

Каковы выводы?

Более крутой тест – тест Колмогорова-Смирнова. Он проверяет на соответствие не нормальному, а любому распределению: \(H_0 \sim F(x)\).

ks.test(w, "pnorm", mean = mean(w), sd = sd(w))
## Warning in ks.test(w, "pnorm", mean = mean(w), sd = sd(w)): ties should not
## be present for the Kolmogorov-Smirnov test
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  w
## D = 0.092203, p-value = 0.5821
## alternative hypothesis: two-sided

Что делать, если вывод функции не понятен? Правильно, ?ks.test.

Теперь мы подошли к третьему способу проверки соответствия распределений: график квантиль-квантиль. Чтобы было проще работать, стандартизуем наши данные (сохраним в новой переменной):

w.st <- (w - mean(w))/sd(w)

И строим наш график. А вот как его читать я расскажу на занятии, ибо тут писать долго и сейчас лень :)

qqnorm(w.st) # Рисуем график
qqline(w.st) # Добавляем линию квартилей

И какие выводы вы можете сделать после этого объяснения?

Резко погружаемся в финансовый мир

В основе бесчисленного числа современных экономических теорий лежит гипотеза о том, что движение котировок подчиняются логнормальному процессу, то есть приращения распределены нормально. На этом строится вся наука о ценообразовании опционов и пресловутый Блэк-Шоулз с их Нобелевской премией.

А давайте проверим эти гипотезы на деле?

Получение финансовых данных

Для нашей тестовой задачи мы возьмем стандартные данные, которые есть в базовом R. Но для настоящего анализа (то есть для ваших домашних работ) вам будут нужны реальные котировки ценных бумаг. Вы вольны брать любые числа, откуда угодно. Но вот две хорошие идеи:

  1. Котировки с сайта Финама.
  2. Котировки иностранных бумаг с сайта Yahoo.finance.

Как красиво получать данные из этих двух источников читайте в самом конце этого файла, а пока давайте анализировать встроенные данные европейских фондовых рынков. Возьмем индекс Британской фондовой биржи FTSE:

prices <- EuStockMarkets[, "FTSE"]
plot(prices, type="l") # type="l" чтобы была линия, а не точки

Теперь нужно перейти от цен к изменениям цен, то есть к приращениям: \[R_t = \frac{P_t-P_{t-1}}{P_t}\] Команда diff вернет нам вектор из приращений в каждой точке. Хотелось бы теперь просто его разделить на значение цен, но нельзя - они разной длины. Поэтому выкинем первый элемент из вектора цен и разделим:

R <- diff(prices)/prices[-1] # Что происходит тут особенно важно понять

Посмотрим на гистограмму и базовые характеристики:

hist(R, breaks = 50)

mean(R)
## [1] 0.000400268
median(R)
## [1] 8.020747e-05
sd(R)
## [1] 0.007951585

Проверим нормальность тестом:

shapiro.test(R)
## 
##  Shapiro-Wilk normality test
## 
## data:  R
## W = 0.98041, p-value = 2.837e-15

И убойным методом:

qqnorm(R)
qqline(R)

Тяжелые хвосты или легкие? Что это значит?

Заключение

Что читать?

Вместо финала дам список базовых книг/статей/курсов по R:

  1. Все основы, очень грамотно и увлекательно, но на английском: Try R course. Очень советую пройти. Готов начислять баллы за бейджики в этом курсе :)
  2. Книга, уже на русском, довольно неплохая: Наглядная статистика. Используем R!
  3. Ну и Google, конечно. Почти все вопросы, которые могут возникнуть у новичков в R уже кто-то где-то обсуждал, так что не стесняйтесь искать.

Что делать? (ДЗ #1)

В качестве первого задания предлагаю получить котировки любой ценной бумаги (или нефти, или чего угодно, что торгуется на рынках) и:

  1. Построить красивую гистограмму доходностей.
  2. Проверить на нормальность двумя статистическими тестами и сделать выводы.
  3. Построить график квантиль-квантиль и сделать выводы.

Формат сдачи – небольшой отчет в свободной форме в любом формате (лучше pdf). Там должно быть Ваше имя, описание того, что Вы делали, Ваши выводы. И обязательно исходный код. Небольшие части кода можно вставлять прямо в отчет, но можно и прислать отдельный R-скрипт. Списывать код – не хорошо. Консультироваться у меня – хорошо.

Примечения

Получение финансовых данных с Yahoo.Finance

Для всего в R есть готовые пакеты, верно? Воспользуемся пакетом quantmod. Чтобы установить пакет, которого у вас нет выполните:

install.packages("quantmod")

После того, как пакет установлен (или если он уже был установлен), то просто выполняем:

library(quantmod) # Эта команда подключает пакет
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## Loading required package: TTR
## Version 0.4-0 included new data defaults. See ?getSymbols.

Теперь у нас появилась пара крутых функций. Нам нужна пока getSymbols – она умеет сама скачивать котировки с сайта Yahoo:

getSymbols("AAPL", from="2015-01-01") # Поулчаем котировки Apple с начала года
##     As of 0.4-0, 'getSymbols' uses env=parent.frame() and
##  auto.assign=TRUE by default.
## 
##  This  behavior  will be  phased out in 0.5-0  when the call  will
##  default to use auto.assign=FALSE. getOption("getSymbols.env") and 
##  getOptions("getSymbols.auto.assign") are now checked for alternate defaults
## 
##  This message is shown once per session and may be disabled by setting 
##  options("getSymbols.warning4.0"=FALSE). See ?getSymbols for more details.
## [1] "AAPL"
prices <- Ad(AAPL) # Функция Ad запирает Adjusted цены
plot(prices)

Поулчение данных с сайта Финам

Идем на сайт экспорта котировок, выбираем нужную бумагу и ставим параметры. Мы пока будем использовать котировки закрытия, то есть одна цена – один торговый день. Важный момент: разделителем лучше ставтье, а формат файла - csv. Скачиваем файл и кладем его в рабочую директорию R (узнать ее поможет getwd() ). Теперь открываем его (я скачал акции АО Газпром):

finam <- read.csv(file = "GAZP.csv", header = TRUE, sep = ";") # Только поменяйте путь до файла.
prices <- finam$X.CLOSE.
plot(prices, type = "l")