Исследование причин повышенного роста цен на первичном рынке жилья в субъектах ДФО
Гришин Данил
Постановка исследовательского вопроса
Начиная с 1 квартала 2018, динамика цен на первичное жилье в ДФО превышает общероссийские показатели (рис. 1). Разрыв в ценах увеличивался вплоть до 3 квартала 2020.
Целью исследовательского проекта является изучение причин повышенного роста цен на первичное жилье в субъектах Дальневосточного федерального округа.
Изучив данные на ресурсах fedstat.ru и cbr.ru, мной был выдвинут следующий исследовательский вопрос: существует ли взаимосвязь динамики цен на первичное жилье в субъектах ДФО от динамики цен на стройматериалы, объема введенного в эксплуатацию жилья, среднесписочной численности рабочих по ВЭД строительство, количества выданных ипотечных жилищных кредитов, средневзвешенной ставки по ипотечному кредитованию и от “Южности” субъекта ДФО.
Согласно данным рисунка 2 сложно выдвинуть гипотезу о влиянии вышеуказанных переменных на динамику цен.
Тем не менее можно предположить следующие исследовательские гипотезы:
- Чем выше динамика цен на стройматериалы в регионе, тем выше динамика цен на жилье;
- Чем выше динамика ввода жилья в эксплуатацию, тем ниже динамика цен на жилье;
- Чем выше динамика объема среднесписочной численности рабочих, тем ниже динамика цен на жилье;
- Чем выше динамика количества выданных ипотечных кредитов, тем выше динамика цен на жилье;
- В южных регионах ДФО динамика цен на жилье в среднем выше, чем в северных регионах;
Наиболее интересна гипотеза, касающаяся влияния средневзвешенной ставки на динамику цен. Как известно, в конце 2019 для дальневосточных регионов была введена программа льготной ипотеки (Дальневосточная ипотека), благодаря которой стало возможным получить ипотеку по 2% годовых. На рис. 3 можно видеть эффект от введения дальневосточной ипотеки. Можно предположить, что, чем выше средневзвешенная ставка, тем ниже динамика цен на жилье.
Как было сказано ранее, для проведения исследования использовались данные ресурсов fedstat.ru и cbr.ru. Для проведения исследования использовались панельные данные по 10 субъектам ДФО за период с начала 2018 по сентябрь 2021.
Используемые переменные и их обозначения представлены в таблице ниже:
| variable | label | Source |
|---|---|---|
| ind_price | Темп роста цен на первичном рынке, квартал к соответствующему кварталу предыдущего года | www.fedstat.ru |
| ind_icp | Индекс цен на продукцию (затраты, услуги) инвестиционного назначения, месяц к соответствующему месяцу предыдущего года | www.fedstat.ru |
| ind_work | Индекс объема введенного жилья, месяц к соответствующему месяцу предыдущего года | www.fedstat.ru |
| ind_emp | Индекс среднесписочной численности рабочих по ВЭД строительство, месяц к соответствующему месяцу предыдущего года | www.fedstat.ru |
| ind_mor | Индекс количества выданных ипотечных жилищных кредитов, месяц к соответствующему месяцу предыдущего года | www.cbr.ru |
| rate | Средневзвешенная ставка по ипотечным жилищным кредитам, предоставленным физическим лицам-резидентам в рублях за месяц, % | www.cbr.ru |
| South | Дамми переменная для обозначения южных регионов ДФО | - |
Следует обратить внимание на то, что зависимая переменная ind_price является квартальной в то время, как остальные факторы - месячные. Для проведения расчетов все месячные переменные были переведены в квартальный формат путем вычисления среднего арифметического (в R с использованием функции aggregate()). Более подробно, как проводилась трасформация данных, можно посмотреть ниже в блоке с кодом.
Итоговую структуру фрейма данных, используемого для вычислений, можно посмотреть далее:
| Date | Subject | ind_price | ind_icp | ind_work | ind_mor | rate | ind_emp | South |
|---|---|---|---|---|---|---|---|---|
| 2018 Q1 | Amur | 100.99 | 103.30333 | 66.46100 | 223.9772 | 9.726667 | 114.16366 | 1 |
| 2018 Q1 | Bur | 101.00 | 107.59667 | 351.89126 | 196.3863 | 9.683333 | 112.89201 | 1 |
| 2018 Q1 | Evr | 91.03 | 100.04333 | 209.20228 | 230.0180 | 9.640000 | 89.18089 | 1 |
| 2018 Q1 | Hab | 99.76 | 102.70000 | 218.99843 | 174.3462 | 9.730000 | 107.60215 | 1 |
| 2018 Q1 | Kam | 99.16 | 123.84333 | 125.27660 | 167.1567 | 9.740000 | 84.67157 | 0 |
| 2018 Q1 | Mag | 95.58 | NA | 150.00000 | 159.8293 | 9.666667 | 112.75920 | 0 |
| 2018 Q1 | Prim | 109.73 | 103.51333 | 150.66497 | 155.2942 | 9.816667 | 100.37176 | 1 |
| 2018 Q1 | Sah | 106.66 | 97.17333 | 69.47011 | 180.0037 | 9.703333 | 122.83973 | 1 |
| 2018 Q1 | Yakut | 101.75 | 100.40333 | 164.52810 | 146.6402 | 9.856667 | 138.86759 | 0 |
| 2018 Q1 | Zab | 100.66 | 102.81667 | 90.77176 | 186.9600 | 9.640000 | 123.36863 | 1 |
Описание данных
Взглянем на диаграммы рассеяния основных факторов (рис. 4 и 5).
## Warning: Removed 7 rows containing non-finite values (stat_smooth).## Warning: Removed 7 rows containing missing values (geom_point).## Warning: Removed 3 rows containing non-finite values (stat_smooth).## Warning: Removed 3 rows containing missing values (geom_point).## Warning: Removed 3 rows containing non-finite values (stat_smooth).## Warning: Removed 3 rows containing missing values (geom_point).## Warning: Removed 3 rows containing non-finite values (stat_smooth).## Warning: Removed 3 rows containing missing values (geom_point).
## Warning: Removed 3 rows containing non-finite values (stat_smooth).## Warning: Removed 3 rows containing missing values (geom_point).
На представленных диаграммах довольно трудно выявить имеющиеся взаимосвязи. Если исключить аутлайеры, то практически во всех случаях будет наблюдаться отсуствие взаимосвязи с динамикой цен на жилье. Можно наблюдать только слабую отрицательную зависимость динамики цен на жилье от количества ипотечных кредитов. А также видно, что с повышением ставки (рис. 5) облако точек на диаграмме становится немного ниже (отрицательная взаимосвязь). Также облако становится плотнее, что говорит о наличии гетероскедастичности.
На рис. 6 также видно, что регионов с более высокой динамикой цен на жилье больше среди южных.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.## Warning: Removed 3 rows containing non-finite values (stat_density).
## Warning: Couldn't find skimmers for class: yearqtr; No user-defined `sfl`
## provided. Falling back to `character`.Взглянем на описательные статистики. В среднем за период с 2018 по 3 квартал 2021 динамика цен на жилье составляла 110,38 % г/г. Стройматериалы дорожали в среднем с темпом 105,75 % г/г. Средняя динамика ввода жилья 135,34 % г/г. Динамика количества ипотечных кредитов 125,03 % г/г. Средняя ставка равна 8,464 % и динамика среднесписочной численности рабочих равна 103,62 % г/г. Выборка состоит из 150 наблюдений, из которых 45 относятся с северным регионам.
| Date | Subject | ind_price | ind_icp | ind_work | ind_mor | rate | ind_emp | South | |
|---|---|---|---|---|---|---|---|---|---|
| Min. :2018 | Bur :15 | Min. : 90.72 | Min. : 92.67 | Min. : 23.33 | Min. : 75.22 | Min. : 5.813 | Min. : 66.58 | 0: 45 | |
| 1st Qu.:2019 | Zab :15 | 1st Qu.:102.19 | 1st Qu.:103.67 | 1st Qu.: 79.74 | 1st Qu.: 98.82 | 1st Qu.: 6.967 | 1st Qu.: 92.96 | 1:105 | |
| Median :2020 | Yakut :15 | Median :109.08 | Median :105.77 | Median :110.43 | Median :121.69 | Median : 9.280 | Median :101.81 | NA | |
| Mean :2020 | Kam :15 | Mean :110.38 | Mean :105.75 | Mean :135.34 | Mean :125.03 | Mean : 8.464 | Mean :103.62 | NA | |
| 3rd Qu.:2021 | Prim :15 | 3rd Qu.:115.79 | 3rd Qu.:107.40 | 3rd Qu.:162.31 | 3rd Qu.:151.30 | 3rd Qu.: 9.662 | 3rd Qu.:112.70 | NA | |
| Max. :2022 | Hab :15 | Max. :197.87 | Max. :131.93 | Max. :509.56 | Max. :230.02 | Max. :10.680 | Max. :164.82 | NA | |
| NA | (Other):60 | NA’s :3 | NA’s :4 | NA | NA | NA | NA | NA |
| Name | df2_no_R_D_South_fact |
| Number of rows | 150 |
| Number of columns | 9 |
| _______________________ | |
| Column type frequency: | |
| character | 1 |
| factor | 2 |
| numeric | 6 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| Date | 0 | 1 | 4 | 7 | 0 | 15 | 0 |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| Subject | 0 | 1 | FALSE | 10 | Bur: 15, Zab: 15, Yak: 15, Kam: 15 |
| South | 0 | 1 | FALSE | 2 | 1: 105, 0: 45 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| ind_price | 3 | 0.98 | 110.38 | 12.54 | 90.72 | 102.19 | 109.08 | 115.79 | 197.87 | ▇▃▁▁▁ |
| ind_icp | 4 | 0.97 | 105.75 | 4.41 | 92.67 | 103.67 | 105.77 | 107.40 | 131.93 | ▁▇▂▁▁ |
| ind_work | 0 | 1.00 | 135.34 | 84.77 | 23.33 | 79.74 | 110.43 | 162.31 | 509.56 | ▇▅▁▁▁ |
| ind_mor | 0 | 1.00 | 125.03 | 31.91 | 75.22 | 98.82 | 121.69 | 151.30 | 230.02 | ▇▇▆▂▁ |
| rate | 0 | 1.00 | 8.46 | 1.51 | 5.81 | 6.97 | 9.28 | 9.66 | 10.68 | ▅▆▁▇▆ |
| ind_emp | 0 | 1.00 | 103.62 | 16.52 | 66.58 | 92.96 | 101.81 | 112.70 | 164.82 | ▂▇▅▂▁ |
Также построим корреляционную матрицу. Между переменными не наблюдается сильной корреляции. Мультиколлинеарность отсутствует.
| ind_price | ind_icp | ind_work | ind_mor | rate | ind_emp | South | |
|---|---|---|---|---|---|---|---|
| ind_price | 1.0000000 | -0.0251390 | 0.3007406 | -0.1635106 | -0.2117199 | 0.0313226 | 0.0725550 |
| ind_icp | -0.0251390 | 1.0000000 | 0.0025257 | -0.0183914 | 0.0776440 | 0.0281000 | -0.1612012 |
| ind_work | 0.3007406 | 0.0025257 | 1.0000000 | -0.1258801 | 0.0375647 | -0.1286240 | -0.0639963 |
| ind_mor | -0.1635106 | -0.0183914 | -0.1258801 | 1.0000000 | -0.2978037 | 0.0983207 | 0.0700318 |
| rate | -0.2117199 | 0.0776440 | 0.0375647 | -0.2978037 | 1.0000000 | 0.1166038 | -0.0042078 |
| ind_emp | 0.0313226 | 0.0281000 | -0.1286240 | 0.0983207 | 0.1166038 | 1.0000000 | -0.0410419 |
| South | 0.0725550 | -0.1612012 | -0.0639963 | 0.0700318 | -0.0042078 | -0.0410419 | 1.0000000 |
Изобразим наши наблюдения в осях первых двух главных компонент
## Warning: ggrepel: 108 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
Оценивание моделей
Оценим простую модель и более сложную:
\[ indprice_i = \beta_1 + \beta_2 indicp_i + \beta_3 indemp_i + u_i \]
\[ indprice_i = \beta_1 + \beta_2 indicp_i + \beta_3 indwork_i + \beta_4 indmor_i + \beta_5 rate_i + \beta_6 indemp_i + \beta_7 South_i + u_i \]
Сравним результаты оценивания моделей:
| (1) | (2) | |
|---|---|---|
| (Intercept) | 116.031 *** | 122.644 *** |
| (25.830) | (24.675) | |
| ind_icp | -0.073 | 0.024 |
| (0.238) | (0.218) | |
| ind_emp | 0.024 | 0.100 |
| (0.063) | (0.058) | |
| ind_work | 0.046 *** | |
| (0.012) | ||
| ind_mor | -0.092 ** | |
| (0.032) | ||
| rate | -2.553 *** | |
| (0.668) | ||
| South1 | 3.158 | |
| (2.144) | ||
| N | 143 | 143 |
| R2 | 0.002 | 0.209 |
| logLik | -563.672 | -547.000 |
| AIC | 1135.343 | 1110.000 |
| *** p < 0.001; ** p < 0.01; * p < 0.05. | ||
Проинтепретируем коэффициенты для сложной модели:
- При увеличении динамики роста цен на стройматериалы на 1 п.п. цены на жилье увеличатся на 0,024 п.п.;
- При увеличении динамики среднесписочной численности рабочих на 1 п.п. цены на жилье увеличатся на 0,1 п.п.;
- При увеличении динамики ввода жилья на 1 п.п., цены на жилье возрастут на 0,046 п.п.;
- При увеличении динамики выданных ипотечных кредитов на 1 п.п. цены на жилье снизятся на 0,092 п.п.;
- Увеличение ставки по ипотеке на 1 п.п. приводит к снижению цен на жилье на 2,553 п.п.;
- При прочих равных условиях динамика цен на жилье выше в южных регионах;
Доверительные интервалы для коэффициентов сложной модели (использованы робастные к гетероскедастичности методы):
## Warning: пакет 'estimatr' был собран под R версии 4.1.2| 2.5 % | 97.5 % | |
|---|---|---|
| (Intercept) | 69.3989143 | 175.8881559 |
| ind_icp | -0.4226827 | 0.4702564 |
| ind_work | -0.0209953 | 0.1125165 |
| ind_mor | -0.1412211 | -0.0436870 |
| rate | -3.6033002 | -1.5018767 |
| ind_emp | 0.0221575 | 0.1782926 |
| South1 | -1.8078917 | 8.1239057 |
Проведем тест Уальда для определения лучшей модели (использованы робастные к гетероскедастичности методы).
| Res.Df | Df | Chisq | Pr(>Chisq) |
|---|---|---|---|
| 140 | |||
| 136 | 4 | 40.4 | 3.62e-08 |
Так как в проведенном тесте p-value < 0.01, то H0 о том, что предпочитается простая модель, отвергается. Следовательно, в нашем исследование предпочтение отдается более сложной модели.
Также проведем тест Рамсея для выбранной модели.
##
## RESET test
##
## data: model2_hc3
## RESET = 0.49751, df1 = 2, df2 = 134, p-value = 0.6092Так как p-value > 0.01, то модель специфицирована верно и пропущенных переменных нет.
Также для сложной модели проведем тест Бройша-Пагана.
##
## studentized Breusch-Pagan test
##
## data: model2
## BP = 27.306, df = 6, p-value = 0.0001269P-value < 0.01, следовательно нулевая гипотеза о том, что квадраты остатков не предсказываются регрессорами, отвергается. В данных есть гетероскедастичность.
Также рассчитаем коэффициенты вздутия дисперсии с помощью функции vif() библиотеки car
## ind_icp ind_work ind_mor rate ind_emp South
## 1.033431 1.034837 1.137570 1.131706 1.052774 1.038991Все коэффициенты находятся вблизи единицы. Поэтому можно предположить отсутствие мультиколлинеарности.
Оценивание LASSO регресии. Сравнение МНК и LASSO по прогнозной силе
Разделим выборку. Для обучения будем использовать 90% имеющихся данные. На оставшихся 10-ти процентах будем проводить тесты.
Посмотрим на коэффициенты оцененных моделей (простая и сложная МНК):
| (1) | (2) | |
|---|---|---|
| (Intercept) | 110.660 *** | 125.330 *** |
| (27.581) | (26.400) | |
| ind_icp | 0.005 | 0.044 |
| (0.255) | (0.229) | |
| ind_emp | -0.005 | 0.112 |
| (0.066) | (0.062) | |
| ind_work | 0.056 *** | |
| (0.013) | ||
| ind_mor | -0.101 ** | |
| (0.034) | ||
| rate | -3.139 *** | |
| (0.745) | ||
| South1 | 2.131 | |
| (2.331) | ||
| N | 128 | 128 |
| R2 | 0.000 | 0.244 |
| logLik | -507.335 | -489.432 |
| AIC | 1022.670 | 994.864 |
| *** p < 0.001; ** p < 0.01; * p < 0.05. | ||
Также посмотрим на коэффициенты при наилучшей LASSO регресии
## 8 x 1 sparse Matrix of class "dgCMatrix"
## s0
## (Intercept) 125.68351481
## (Intercept) .
## ind_icp 0.04032434
## ind_work 0.05606855
## ind_mor -0.10017198
## rate -3.12133610
## ind_emp 0.11016155
## South1 2.08243328Теперь построим прогноз на тестовую выборку по каждой модели и сравним качество полученных прогнозов с помощью RSS:
| y fact | y pred ols sh | y pred ols lg | y pred LASSO | |
|---|---|---|---|---|
| 135 | 102.34 | 110.6245 | 117.6673 | 117.6767 |
| 136 | 100.00 | 110.7096 | 105.8074 | 105.8787 |
| 137 | 112.08 | 110.6557 | 114.0889 | 114.0431 |
| 138 | 104.83 | 110.5532 | 120.2168 | 120.1638 |
| 139 | 111.63 | 110.6341 | 109.9351 | 109.9822 |
| 140 | 116.45 | 110.5311 | 113.6376 | 113.5983 |
| 141 | 122.39 | 110.4861 | 135.4644 | 135.3075 |
| 142 | 130.55 | 110.5450 | 118.4602 | 118.4380 |
| 144 | 116.06 | 110.6825 | 113.9118 | 113.8907 |
| 145 | 100.21 | 110.6400 | 111.5469 | 111.5610 |
| 146 | 100.07 | 110.6388 | 128.2039 | 128.1495 |
| 147 | 108.12 | 110.6541 | 122.0652 | 121.9874 |
| 148 | 109.32 | 110.5080 | 121.3848 | 121.3100 |
| 149 | 111.48 | 110.6328 | 115.3575 | 115.3649 |
| 150 | 131.51 | 110.4867 | 122.2804 | 122.1826 |
| RSS_ols_sh | 1495.973 |
| RSS_ols_lg | 2202.234 |
| RSS_lasso_best | 2193.330 |
Судя по RSS, наилучшей моделью оказалась простая МНК регрессия.
Построение прогноза для финальной модели
Созданим новое наблюдение и построим прогноз для него, а также предиктивный интервал.
| ind_icp | ind_emp |
|---|---|
| 107.9 | 105.9 |
Как можно видеть новое наблюдение было создано для региона, у которого:
- Динамика цен на стройматериалы составляет 107,9 % г/г;
- Динамика среднесписочной численности рабочих по ВЭД строительство 105,9 % г/г;
Для такого индивида прогнозное значение его дохода составляет 110.61 рублей. А доверительный интервал равен
| fit | lwr | upr |
|---|---|---|
| 110.6099 | 84.97789 | 136.242 |
В целом прогноз не очень точный, так как очень велика доля неопределенности (цены на жилье все-таки пойдут вверх или вниз?).
Блок кода
# Загрузка библиотек
library(readxl)
library(tidyverse) # манипуляции с данными
library(skimr) # описательные статистики
library(tidymodels) # удобства для оценивания моделей
library(sjPlot) # визуализация коэффициентов
library(car) # проверка сложных гипотез
library(lmtest) # тесты для моделей
library(texreg) # таблички для сравнения моделей
library(rio) # import / export всех форматов
library(labelled)
library(sjlabelled)
library(ggplot2)
library(plotly)
library(dplyr)
library(huxtable) # несколько таблиц вместе
library(corrplot) # картинки для корреляционной матрицы
library(tseries)
library(reshape)
library(knitr)
# Выгружаем датасет
ind_price <- read_excel("data.xlsx", sheet = "ind_price") # y
ind_icp <- read_excel("data.xlsx", sheet = "ind_icp") #x1
ind_work <- read_excel("data.xlsx", sheet = "ind_work") #x2
ind_emp <- read_excel("data.xlsx", sheet = "ind_emp") #x3
ind_mor <- read_excel("data.xlsx", sheet = "ind_mor") #x4
rate <- read_excel("data.xlsx", sheet = "rate") #x5
# Выгружаем обозначения переменных
Labels <- read_excel("data.xlsx", sheet = "Labels")
Labels <- rbind(Labels, data.frame(variable = 'South', label = 'Дамми переменная для обозначения южных регионов ДФО', Source = '-'))
kable(Labels, escape=FALSE)
# kable(glimpse(df2[0:10,]), escape=FALSE)
# Создадим вектор субъектов для замены названий столбцов
regs <- c('Russia','DFO','Bur','Zab','Yakut','Kam','Prim','Hab','Amur','Mag','Sah','Evr')
### ind_price - y
y <- as.data.frame(t(ind_price[-1]))
colnames(y) <- regs
qtrs <- ts(y, start = c(2010,1), frequency = 4)
qtrss <- as.yearqtr(time(qtrs), format = "%Y-%q")
y['Qtr'] <- qtrss
y_resh <- melt(y, id=c("Qtr"))
colnames(y_resh) <- c('Date','Subject','ind_price')
### ind_icp - x1
x1 <- as.data.frame(t(ind_icp[-1]))
colnames(x1) <- regs
mons <- ts(x1, start = c(2010,1), frequency = 12)
monss <- as.yearmon(time(mons), format = "%Y-%m")
x1['Month'] <- monss
qu <- quarters(x1$Month)
ye <- strftime(x1$Month, format = "%Y")
x1['qu'] <- qu
x1['ye'] <- ye
x1_qtr <- aggregate(cbind(x1[1:12]),list(x1$qu,x1$ye), mean)
x1_qtr <- x1_qtr[,-c(1,2)]
x1_qtr['Qtr'] <- qtrss
x1_qtr_resh <- melt(x1_qtr, id=c("Qtr"))
colnames(x1_qtr_resh) <- c('Date','Subject','ind_icp')
### ind_work - x2
x2 <- as.data.frame(t(ind_work[-1]))
colnames(x2) <- regs
mons <- ts(x2, start = c(2010,1), frequency = 12)
monss <- as.yearmon(time(mons), format = "%Y-%m")
x2['Month'] <- monss
qu <- quarters(x2$Month)
ye <- strftime(x2$Month, format = "%Y")
x2['qu'] <- qu
x2['ye'] <- ye
x2_qtr <- aggregate(cbind(x2[1:12]),list(x2$qu,x2$ye), mean)
x2_qtr <- x2_qtr[,-c(1,2)]
x2_qtr['Qtr'] <- qtrss
x2_qtr_resh <- melt(x2_qtr, id=c("Qtr"))
colnames(x2_qtr_resh) <- c('Date','Subject','ind_work')
### ind_emp - x3
x3 <- as.data.frame(t(ind_emp[-1]))
colnames(x3) <- regs
mons <- ts(x3, start = c(2018,1), frequency = 12)
monss <- as.yearmon(time(mons), format = "%Y-%m")
x3['Month'] <- monss
qu <- quarters(x3$Month)
ye <- strftime(x3$Month, format = "%Y")
x3['qu'] <- qu
x3['ye'] <- ye
x3_qtr <- aggregate(cbind(x3[1:12]),list(x3$qu,x3$ye), mean)
x3_qtr <- x3_qtr[,-c(1,2)]
quart <- ts(x3_qtr, start = c(2018,1), frequency = 4)
quarts <- as.yearqtr(time(quart), format = "%Y-%q")
x3_qtr['Qtr'] <- quarts
x3_qtr_resh <- melt(x3_qtr, id=c("Qtr"))
colnames(x3_qtr_resh) <- c('Date','Subject','ind_emp')
### ind_mor - x4
x4 <- as.data.frame(t(ind_mor[-1]))
colnames(x4) <- regs
mons <- ts(x4, start = c(2010,1), frequency = 12)
monss <- as.yearmon(time(mons), format = "%Y-%m")
x4['Month'] <- monss
qu <- quarters(x4$Month)
ye <- strftime(x4$Month, format = "%Y")
x4['qu'] <- qu
x4['ye'] <- ye
x4_qtr <- aggregate(cbind(x4[1:12]),list(x4$qu,x4$ye), mean)
x4_qtr <- x4_qtr[,-c(1,2)]
x4_qtr['Qtr'] <- qtrss
x4_qtr_resh <- melt(x4_qtr, id=c("Qtr"))
colnames(x4_qtr_resh) <- c('Date','Subject','ind_mor')
### rate - x5
x5 <- as.data.frame(t(rate[-1]))
colnames(x5) <- regs
mons <- ts(x5, start = c(2010,1), frequency = 12)
monss <- as.yearmon(time(mons), format = "%Y-%m")
x5['Month'] <- monss
qu <- quarters(x5$Month)
ye <- strftime(x5$Month, format = "%Y")
x5['qu'] <- qu
x5['ye'] <- ye
x5_qtr <- aggregate(cbind(x5[1:12]),list(x5$qu,x5$ye), mean)
x5_qtr <- x5_qtr[,-c(1,2)]
x5_qtr['Qtr'] <- qtrss
x5_qtr_resh <- melt(x5_qtr, id=c("Qtr"))
colnames(x5_qtr_resh) <- c('Date','Subject','rate')
df1 <- cbind(y_resh,x1_qtr_resh[3],x2_qtr_resh[3],x4_qtr_resh[3],x5_qtr_resh[3])
df2 <- merge(df1,x3_qtr_resh,join = 'left',fill='NA')
South_regs <- c('Bur','Zab','Prim','Hab','Amur','Sah','Evr')
df1['South'] <- case_when(df1$Subject == 'Bur' ~ 1,
df1$Subject == 'Zab' ~ 1,
df1$Subject == 'Prim' ~ 1,
df1$Subject == 'Hab' ~ 1,
df1$Subject == 'Amur' ~ 1,
df1$Subject == 'Sah' ~ 1,
df1$Subject == 'Evr' ~ 1)
library("imputeTS")
df1['South'] <- na.replace(df1['South'],0)
df2['South'] <- case_when(df2$Subject == 'Bur' ~ 1,
df2$Subject == 'Zab' ~ 1,
df2$Subject == 'Prim' ~ 1,
df2$Subject == 'Hab' ~ 1,
df2$Subject == 'Amur' ~ 1,
df2$Subject == 'Sah' ~ 1,
df2$Subject == 'Evr' ~ 1)
df2['South'] <- na.replace(df2['South'],0)
# Рис 1 График-сравнение динамики цен в целом по РФ и по ДФО
df2_f <- filter(df2,Subject == c('Russia') | Subject == c('DFO'))
ggplot(df2_f) +
geom_line(mapping = aes(x = Date, y = ind_price, colour = Subject),
size = 1.5) +
theme_classic() +
scale_color_manual(name = "",values = c('#33CCCC','#FF9966'),labels = c("Россия", "ДФО")) +
ylab('') +
xlab('') +
labs(title = "Рис. 1. С 1Q 2018 увеличился разрыв динамики цен на первичное жилье \n в ДФО от общероссийских показателей",
subtitle = "% прироста г/г") +
theme(
plot.title = element_text(face="bold",colour = "black"),
legend.text = element_text(size = 14),
legend.title = element_text(size = 14,face = "bold"),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 14)
)
# Рис 2 - сравнение динамики цен на жилье с ключевыми факторами
df2_f_DFO <- filter(df2, Subject == 'DFO')
df2_f_DFO <- melt(df2_f_DFO, id=c("Date","Subject"))
df_graph <- filter(df2_f_DFO, variable == 'ind_price'
| variable == 'ind_icp'
| variable == 'ind_work'
| variable == 'ind_emp'
| variable == 'ind_mor')
ggplot(df_graph) +
geom_line(mapping = aes(x = Date, y = value, colour = variable), size = 1.5) +
theme_classic() +
scale_color_manual(name = "",values = c('#FF9966','forestgreen','steelblue','#9999FF','orange'),labels = c("Цены на жилье",
"Цены на стройматериалы",
"Объем ввода жилья",
"Ипотечные кредиты",
"Рабочая сила")) +
ylab('') +
xlab('') +
labs(title = "Рис. 2. Динамика ключевых факторов, определяющих динамику цены \n на первичное жилье в ДФО",
subtitle = "% прироста г/г") +
theme(
plot.title = element_text(face="bold",colour = "black"),
legend.text = element_text(size = 12),
legend.title = element_text(size = 12,face = "bold"),
axis.title.x = element_text(size = 12),
axis.title.y = element_text(size = 12),
plot.margin = margin(2, -125, 2, 2)
)
# Рис. 3 динамика ставка по ипотеке (сравнение России и ДФО)
df2_f2 <- filter(df2, Subject == 'Russia' | Subject == 'DFO')
df2_f2_t <- melt(df2_f2, id=c("Date","Subject"))
df_graph2 <- filter(df2_f2_t, variable == 'rate')
ggplot(df_graph2) +
geom_line(mapping = aes(x = Date, y = value, colour = Subject),size = 1.5)+
theme_classic() +
scale_color_manual(name = "",values = c('#33CCCC','#FF9966'),labels = c("Россия", "ДФО")) +
ylab('') +
xlab('') +
labs(title = "С 1Q 2020 ставка по ипотеке в ДФО значительно снизилась \n по сравнению с общероссийскими показателями",
subtitle = "Средневзвешенная ставка, %")+
theme(
plot.title = element_text(face="bold",colour = "black"),
legend.text = element_text(size = 14),
legend.title = element_text(size = 14,face = "bold"),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 14)
)
# Рис 4
library(gridExtra)
df2_South_fact <- mutate_at(df2, vars(South), factor)
df2_clear <- filter(df2_South_fact,Subject != 'Russia', Subject != 'DFO')
# Стройматериалы
g1 <- ggplot(df2_clear) +
geom_point(mapping = aes(x = ind_icp, y = ind_price)) +
stat_smooth(mapping = aes(x = ind_icp, y = ind_price), se = FALSE, colour = 'red',
linetype = 'dotted',size = 1.2) +
labs(title = "Рис. 4") +
theme_classic() +
xlab('Цены на стройматериалы') +
ylab('Цены на жилье') +
labs(subtitle = "% прироста г/г")
# Объем ввода
g2 <- ggplot(df2_clear) +
geom_point(mapping = aes(x = ind_work, y = ind_price)) +
stat_smooth(mapping = aes(x = ind_work, y = ind_price), se = FALSE, colour = 'red',
linetype = 'dotted',size = 1.2) +
theme_classic() +
xlab('Объем ввода жилья') +
ylab('Цены на жилье') +
labs(subtitle = "% прироста г/г")
# Количество ипотек
g3 <- ggplot(df2_clear) +
geom_point(mapping = aes(x = ind_mor, y = ind_price)) +
stat_smooth(mapping = aes(x = ind_mor, y = ind_price), se = FALSE, colour = 'red',
linetype = 'dotted',size = 1.2) +
theme_classic() +
xlab('Количество ипотечных кредитов') +
ylab('Цены на жилье') +
labs(subtitle = "% прироста г/г")
# Рабочая сила
g4 <- ggplot(df2_clear) +
geom_point(mapping = aes(x = ind_emp, y = ind_price)) +
stat_smooth(mapping = aes(x = ind_emp, y = ind_price), se = FALSE, colour = 'red',
linetype = 'dotted',size = 1.2) +
theme_classic() +
xlab('Среднесписочная численность рабочих') +
ylab('Цены на жилье') +
labs(subtitle = "% прироста г/г")
# Ставка
g5 <- ggplot(df2_clear) +
geom_point(mapping = aes(x = rate, y = ind_price)) +
stat_smooth(mapping = aes(x = rate, y = ind_price), se = FALSE, colour = 'red',
linetype = 'dotted',size = 1.2) +
labs(title = "Рис. 5") +
theme_classic() +
xlab('Средневзвешенная ставка по ипотеке, %') +
ylab('Цены на жилье') +
labs(subtitle = "% прироста г/г")
grid.arrange(g1, g2, g3, g4)
# Рис 5 ставка
g5
# Рис 6 South
qplot(data = df2_clear, x = ind_price, fill = South,
geom = 'density', alpha = 0.3) +
xlab('Цены на жилье') +
ggtitle("") +
labs(title = "Рис. 6") +
theme(
plot.title = element_text(face="bold",colour = "black",),
legend.text = element_text(size = 14),
legend.title = element_text(size = 14,face = "bold"),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 14)
) +
guides(alpha = FALSE) +
theme_classic()
df2_no_R_D <- filter(df2,Subject != "Russia", Subject != "DFO")
df2_no_R_D_na <- na.omit(df2_no_R_D)
corr_mat = cor(df2_no_R_D_na[,c(3:9)])
corr_mat
corrplot.mixed(corr_mat, order = 'hclust')
# Изображение наблюдений в осях первых двух главных компонент
library(factoextra)
df2_no_R_D_no_ind_price <- select(df2_no_R_D,-ind_price)
spisok_reg_date <- c()
for (i in seq_along(df2_no_R_D_no_ind_price$Subject)){
spisok_reg_date[i] <- paste(df2_no_R_D_no_ind_price$Subject[i], df2_no_R_D_no_ind_price$Date[i],collapse = " ")
}
rownames(df2_no_R_D_no_ind_price) <- spisok_reg_date
df2_no_R_D_no_ind_price <- na.omit(select(df2_no_R_D_no_ind_price, -Date, -Subject))
df2_no_R_D_no_ind_price_pca <- prcomp(df2_no_R_D_no_ind_price, scale = TRUE)
fviz_eig(df2_no_R_D_no_ind_price_pca)
fviz_pca_biplot(df2_no_R_D_no_ind_price_pca, repel = TRUE)
df2_no_R_D_South_fact <- mutate_at(df2_no_R_D, vars(South), factor)
# Описательные статистики
df2_summary <- summary(df2_no_R_D_South_fact)
df2_summary
write.csv(df2_summary, 'df2_summary.csv',col.names = TRUE, row.names = TRUE)
# Описательные статистики
df2_skim <- skim(df2_no_R_D_South_fact)
df2_skim
write.csv(df2_skim, 'df2_skim.csv',col.names = TRUE, row.names = TRUE)
model1 <- lm(data = df2_no_R_D_South_fact, ind_price ~ ind_icp + ind_emp)
model2 <- lm(data = df2_no_R_D_South_fact, ind_price ~ ind_icp + ind_work + ind_mor + rate + ind_emp + South)
summary(model1)
summary(model2)
huxreg(model1, model2)
confint(model2)
bptest(model2)
# H0: квадраты остатков не предсказываются регрессорами
# H0 отвергается
library(estimatr)
model1_hc3 <- lm_robust(data = df2_no_R_D_South_fact, ind_price ~ ind_icp + ind_emp, se_type = 'HC3')
model2_hc3 <- lm_robust(data = df2_no_R_D_South_fact, ind_price ~ ind_icp + ind_work + ind_mor +
rate + ind_emp + South, se_type = 'HC3')
summary(model1_hc3)
summary(model2_hc3)
huxreg(model1_hc3, model2_hc3)
confint(model2_hc3)
# model1_hc3 vs model2_hc3, F-test
waldtest(model1_hc3, model2_hc3)
# P-value < alpha (=0.01)
# H0 отвергается. Следовательно, предпочитается более сложная модель
resettest(model2_hc3)
# P-value > 0.01
# Модель специфицирована верно
# Пропущенных переменных нет
vif(model2)
## Сравнение с LASSO
library(glmnet)
df2_no_R_D_South_fact_no_na <- na.omit(df2_no_R_D_South_fact)
train_data <- df2_no_R_D_South_fact_no_na[0:128,] # 90% данных для обучения
test_data <- df2_no_R_D_South_fact_no_na[129:nrow(df2_no_R_D_South_fact_no_na),]
y_train <- train_data$ind_price
X0_train <- model.matrix(data = train_data, ind_price ~ ind_icp + ind_work + ind_mor +
rate + ind_emp + South)
y_test <- test_data$ind_price
X0_test <- model.matrix(data = test_data, ind_price ~ ind_icp + ind_work + ind_mor +
rate + ind_emp + South)
lambda = seq(from = 50, to = 0.01, length = 30)
lambda
m_ols_sh = lm(data = train_data, ind_price ~ ind_icp + ind_emp)
m_ols_lg = lm(data = train_data, ind_price ~ ind_icp + ind_work + ind_mor + rate + ind_emp + South)
m_lasso = glmnet(X0_train, y_train, alpha = 1, lambda = lambda)
lasso.cv.out = cv.glmnet(X0_train,y_train,alpha=1)
lasso.bestlam = lasso.cv.out$lambda.min
m_lasso_best = glmnet(X0_train, y_train, alpha = 1, lambda = lasso.bestlam)
coef(m_lasso_best)
predict_ols_sh <- predict(m_ols_sh,test_data)
predict_ols_lg <- predict(m_ols_lg,test_data)
predict_lasso_best <- predict(m_lasso_best,X0_test)
predicts_df <- cbind(test_data$ind_price,predict_ols_sh,predict_ols_lg,predict_lasso_best)
colnames(predicts_df) <- c('y fact','y pred ols sh','y pred ols lg','y pred LASSO')
kable(predicts_df, escape=FALSE)
RSS_ols_sh <- sum((y_test-predict_ols_sh)^2)
RSS_ols_lg <- sum((y_test-predict_ols_lg)^2)
RSS_lasso_best <- sum((y_test-predict_lasso_best)^2)
RSS_df <- rbind(RSS_ols_sh,RSS_ols_lg,RSS_lasso_best)
kable(RSS_df, escape=FALSE)
# Прогноз для финальной модели m_ols_sh
new = tibble(ind_icp = 107.9, ind_emp = 105.9)
new
prognoz = predict(m_ols_sh, newdata = new)
as.numeric(prognoz)
prognoz_int = predict(m_ols_sh, newdata = new,
interval = 'prediction')
prognoz_int