Просмотр кода страницы Анализ временных рядов/R

Анализ временных рядов и регрессионные методы

# Практические примеры в R: Анализ временных рядов и регрессионные методы

## Полный рабочий пример 1: Segregation Model

### Симуляция данных

```r
# Симулируем данные, похожие на реальные результаты Segregation Model
# Параметр: %-SIMILAR-WANTED = 30%

set.seed(123)

# Генерируем нестационарный ряд (с трендом и шумом)
time_steps <- 1:100
true_trend <- 50 + 0.3 * time_steps  # тренд: медленный рост
# добавляем шум и эффект насыщения
noise <- rnorm(100, mean = 0, sd = 2)
saturation <- 100 * (1 - exp(-0.05 * time_steps))  # асимптотический рост
happy_agents_pct <- pmin(true_trend + noise + saturation - 50, 95)

# Визуализируем
plot(time_steps, happy_agents_pct, type='l', 
     main='Segregation Model: % счастливых агентов',
     xlab='Шаг времени', ylab='% счастливых',
     lwd=2, col='darkblue')
grid()
```

### Шаг 1: Проверка на стационарность

```r
# Используем ADF тест (Augmented Dickey-Fuller)
library(tseries)

# Тест на исходном ряде
adf_original <- adf.test(happy_agents_pct)
cat("ADF тест на исходном ряде:\n")
print(adf_original)

# Если p-value > 0.05, ряд нестационарен
# Обычно p-value < 0.05 для стационарного ряда (отвергаем H0 о нестационарности)
```

### Шаг 2: Линейная регрессия и детерендизация

```r
# Строим линейную регрессию для извлечения тренда
regression_model <- lm(happy_agents_pct ~ time_steps)
summary(regression_model)

# Вывод покажет:
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)    
# (Intercept) 50.1234    0.5678  88.27   <2e-16 ***
# time_steps   0.2987    0.0098  30.49   <2e-16 ***

# Извлекаем остатки (детерендированный ряд)
residuals_detrended <- residuals(regression_model)

# Визуализируем остатки
par(mfrow=c(2,2))

plot(time_steps, residuals_detrended, type='l',
     main='Остатки после удаления тренда',
     xlab='Шаг времени', ylab='Остатки')
abline(h=0, lty=2, col='red')

# Гистограмма остатков
hist(residuals_detrended, breaks=20, 
     main='Распределение остатков',
     xlab='Значение', freq=FALSE)
curve(dnorm(x, mean=mean(residuals_detrended), 
            sd=sd(residuals_detrended)), 
      add=TRUE, col='red', lwd=2)

# Q-Q plot
qqnorm(residuals_detrended, main='Q-Q диаграмма')
qqline(residuals_detrended, col='red')

par(mfrow=c(1,1))
```

### Шаг 3: Критерий Дарбина-Уотсона

```r
library(lmtest)

# Основной тест
dw_result <- dwtest(regression_model)
cat("\n=== Тест Дарбина-Уотсона ===\n")
print(dw_result)

# Ручной расчёт для проверки
n <- length(residuals_detrended)
dw_manual <- sum(diff(residuals_detrended)^2) / sum(residuals_detrended^2)
cat("DW статистика (ручной расчёт):", dw_manual, "\n")

# Интерпретация:
# Если DW ~ 2: нет автокорреляции ✓
# Если DW ~ 0: положительная автокорреляция (остатки похожи)
# Если DW ~ 4: отрицательная автокорреляция (остатки чередуются)

if (dw_result$p.value > 0.05) {
  cat("Вывод: Нет значимой автокорреляции (p > 0.05) ✓\n")
} else {
  cat("Вывод: Есть значимая автокорреляция (p < 0.05) ⚠\n")
}
```

### Шаг 4: Автокорреляционная функция (ACF)

```r
# Визуализация ACF для остатков
par(mfrow=c(2,1))

acf(residuals_detrended, main='ACF остатков (исходные)')
pacf(residuals_detrended, main='Частная ACF остатков')

par(mfrow=c(1,1))

# Интерпретация:
# - Если столбцы внутри пунктирных границ → нет значимой автокорреляции
# - Если выходят за границы → есть автокорреляция на данном лаге
```

### Шаг 5: Альтернативный метод - Дифференцирование

```r
# Первые разности (простое дифференцирование)
diff_happy <- diff(happy_agents_pct)

# Визуализируем
par(mfrow=c(2,1))

plot(happy_agents_pct, type='l', main='Исходный ряд')
plot(diff_happy, type='l', main='Первые разности')

par(mfrow=c(1,1))

# Проверяем стационарность разностей
adf_diff <- adf.test(diff_happy)
cat("ADF тест на дифференцированном ряде:\n")
print(adf_diff)

# Если p-value < 0.05, ряд стационарен ✓
```

---

## Полный рабочий пример 2: Minority Game

### Специфика: высокая волатильность

```r
# Minority Game: процент успешных выборов агента во времени
# Данные более шумные, чем в Segregation

set.seed(456)
time_steps_mg <- 1:150

# Создаём более волатильный ряд
base_success <- 50 + 2 * sin(time_steps_mg / 20)  # слабый циклический тренд
noise_mg <- rnorm(150, mean=0, sd=4)
success_rate <- pmax(pmin(base_success + noise_mg, 90), 10)  # bounds: 10-90%

plot(time_steps_mg, success_rate, type='l',
     main='Minority Game: % успешных выборов',
     xlab='Раунд игры', ylab='Успешность (%)',
     lwd=2, col='darkgreen')
grid()
```

### Анализ: Проверка стационарности

```r
# ADF тест
adf_mg <- adf.test(success_rate)
cat("ADF тест:\n")
print(adf_mg)

# KPSS тест (альтернативный)
library(tseries)
kpss_test <- kpss.test(success_rate)
cat("\nKPSS тест:\n")
print(kpss_test)

# Примечание: если ADF не отвергает H0 (p > 0.05) 
# и KPSS отвергает H0 (p < 0.05), то ряд точно нестационарен
```

### Анализ: Более сложная регрессия

```r
# В Minority Game может быть эффект памяти
# Проверим: зависит ли успех в t от успеха в t-1

# Создаём лаговые переменные
lag1_success <- c(NA, success_rate[-150])
lag2_success <- c(NA, NA, success_rate[-149:150])

# Удаляем NA
valid_idx <- 3:150
y <- success_rate[valid_idx]
x_lag1 <- lag1_success[valid_idx]
x_lag2 <- lag2_success[valid_idx]
time_x <- time_steps_mg[valid_idx]

# Регрессия с лагами
ar_model <- lm(y ~ time_x + x_lag1 + x_lag2)
summary(ar_model)

# Проверяем остатки
residuals_ar <- residuals(ar_model)
dw_ar <- dwtest(ar_model)
print(dw_ar)

# Если лаги значимы, это указывает на автокорреляцию в исходных данных
```

---

## Пример 3: Сравнение методов детерендизации

```r
# Используем данные Segregation Model из примера 1

# Метод 1: Линейная регрессия
residuals_linear <- residuals(lm(happy_agents_pct ~ time_steps))

# Метод 2: Дифференцирование
residuals_diff <- diff(happy_agents_pct)

# Метод 3: Полиномиальная регрессия (более гибкая)
poly_model <- lm(happy_agents_pct ~ poly(time_steps, 2))
residuals_poly <- residuals(poly_model)

# Визуализация результатов
par(mfrow=c(2,2))

plot(time_steps, residuals_linear, type='l', 
     main='Метод 1: Линейная регрессия',
     ylab='Остатки')
abline(h=0, col='red', lty=2)

plot(time_steps[-1], residuals_diff, type='l',
     main='Метод 2: Дифференцирование',
     ylab='Первые разности')
abline(h=0, col='red', lty=2)

plot(time_steps, residuals_poly, type='l',
     main='Метод 3: Полиномиальная регрессия',
     ylab='Остатки')
abline(h=0, col='red', lty=2)

# Сравнение качества
cat("=== Сравнение качества моделей ===\n")
cat("Линейная регрессия:\n")
print(summary(lm(happy_agents_pct ~ time_steps))$adj.r.squared)

cat("Полиномиальная регрессия:\n")
print(summary(poly_model)$adj.r.squared)

par(mfrow=c(1,1))
```

---

## Пример 4: Экспорт данных из NetLogo в R

### Файл NetLogo (.csv экспорт)

```r
# Допустим, вы экспортировали данные из NetLogo как 'segregation_data.csv'
# Структура файла:
# [step] [happy-percent] [num-unhappy]
# 0      50.2            245
# 1      52.1            231
# ...

# Загружаем данные
data_nl <- read.csv('segregation_data.csv')

# Проверяем структуру
head(data_nl)
str(data_nl)

# Вынимаем колонку счастья
happy_ts <- data_nl$happy_percent

# Теперь применяем все методы анализа из примеров выше
model <- lm(happy_ts ~ data_nl$step)
summary(model)

dw <- dwtest(model)
print(dw)
```

---

## Полезные функции и пакеты

| Задача | Функция | Пакет |
|--------|---------|-------|
| ADF тест | `adf.test()` | tseries |
| KPSS тест | `kpss.test()` | tseries |
| Дарбина-Уотсона | `dwtest()` | lmtest |
| ACF/PACF | `acf()`, `pacf()` | stats (встроен) |
| Дифференцирование | `diff()` | stats |
| Линейная регрессия | `lm()` | stats |
| Полиномиальная | `lm(y ~ poly(x,2))` | stats |
| Визуализация | `plot()`, `lines()` | graphics |
| Сглаживание (MA) | `filter()` | stats |

---

## Интерпретация результатов: Таблица принятия решений

| Ситуация | DW значение | p-value | Действие |
|----------|-------------|---------|----------|
| Идеально | 1.8-2.2 | > 0.05 | Используйте модель как есть ✓ |
| Положительная автокор. | < 1.5 | < 0.05 | Добавьте лаги переменных |
| Отрицательная автокор. | > 2.5 | < 0.05 | Проверьте спецификацию модели |
| Неопределённость | 1.5-1.8 | ? | Используйте дополнительные тесты (ACF) |

---

## Дополнительные ресурсы

- Документация `lmtest`: `?dwtest`
- Документация `tseries`: `?adf.test`
- Для ARIMA моделей: пакет `forecast`
- Для более сложных: пакет `dynlm` (Dynamic Linear Models)