1. Wprowadzenie i Uzasadnienie Wyboru Zjawiska

Zjawisko badane: Przestrzenne zróżnicowanie cen mieszkań w Polsce.

Uzasadnienie: Ceny nieruchomości mieszkalnych stanowią jedno z kluczowych zjawisk ekonomicznych o wyraźnym wymiarze przestrzennym. Zgodnie z pierwszym prawem geografii Toblera („wszystko jest powiązane ze wszystkim, ale bliskie obiekty są bardziej powiązane ze sobą niż odległe”), można oczekiwać, że rynki mieszkaniowe w geograficznie bliskich miastach wykazują podobne tendencje cenowe – zjawisko przestrzennej autokorelacji.

Hipoteza badawcza: Ceny mieszkań w polskich miastach wykazują istotną dodatnią autokorelację przestrzenną – tj. wysokie ceny w jednym mieście są powiązane z wysokimi cenami w miastach geograficznie bliskich (i analogicznie dla cen niskich).

Dane: Zbiór ofert sprzedaży mieszkań z portalu Otodom Analytics (czerwiec 2024 r.), zawierający ~21 000 ogłoszeń z 15 największych polskich miast.

Źródło danych: [https://www.kaggle.com/datasets/krzysztofjamroz/apartment-prices-in-poland/?select=apartments_pl_2023_08.csv]

1.1 Determinanty Cen Mieszkań

Na dane pozwoliły zidentyfikować kluczowe determinanty cenowe:

Zmienna Spodziewany kierunek wpływu Uzasadnienie
Metraż (squareMeters) (+) całkowita cena, (−) cena/m² Efekt skali przy zakupie dużych lokali
Liczba pokoi (rooms) (+) Większy lokal = większa wartość użytkowa
Odległość od centrum (centreDistance) (−) Klasyczna renta lokalizacyjna
Liczba punktów POI (poiCount) (+) Bogata infrastruktura usługowa podnosi wartość
Obecność windy (hasElevator) (+) Udogodnienie zwiększające dostępność i komfort
Parking (hasParkingSpace) (+) Dodatkowy zasób zwiększający wartość
Rok budowy (buildYear) (+) Nowsze budownictwo = wyższy standard

2. Przygotowanie Danych

2.1 Czyszczenie i Selekcja

Po wczytaniu surowego zbioru (21501 obserwacji) wykonano następujące operacje:

  1. Selekcja zmiennych – wybrano kolumny istotne dla modelu.
  2. Rekodowanie – zmienne binarne (hasElevator, hasParkingSpace) przekonwertowano z formatu yes/no na 1/0.
  3. Usunięcie braków danych – z analizy wykluczono rekordy z wartościami NA.
  4. Filtracja cen – usunięto ogłoszenia z zerową lub ujemną ceną.
  5. Agregacja – zbiór zagregowano do poziomu 15 miast (jedna obserwacja = jedno miasto).

Po czyszczeniu: 17554 obserwacji w 15 miastach.

2.2 Tabela Zagregowanych Statystyk

Zagregowane statystyki dla 15 miast (posortowane wg śr. ceny/m²)
Miasto Śr. cena/m² Śr. m² Śr. pokoje Śr. odl. centrum Śr. POI Śr. winda Śr. parking Śr. rok bud. Szer. geogr. Dług. geogr. N ogłoszeń
warszawa 18517.50 57.78 2.61 5.94 22.92 0.69 0.34 1988.84 52.23 21.01 6140
krakow 17154.68 55.25 2.55 4.19 21.81 0.50 0.33 1990.17 50.07 19.95 2712
gdansk 15450.79 59.01 2.73 4.81 18.15 0.55 0.12 1995.55 54.35 18.61 1706
gdynia 14108.05 63.78 2.93 4.15 12.51 0.46 0.14 1990.02 54.51 18.51 739
wroclaw 13649.09 56.57 2.60 3.87 18.93 0.51 0.43 1988.08 51.11 17.03 1927
poznan 11661.13 57.68 2.63 3.60 19.54 0.51 0.42 1983.01 52.41 16.92 569
rzeszow 10935.30 60.10 2.93 2.37 12.46 0.54 0.36 2003.05 50.04 22.00 127
bialystok 10213.47 51.61 2.58 2.06 18.89 0.38 0.19 1991.98 53.14 23.15 193
lublin 10064.68 58.83 2.86 3.01 17.62 0.40 0.23 1987.95 51.24 22.55 498
szczecin 9702.92 63.26 2.80 3.58 21.28 0.30 0.34 1972.23 53.43 14.54 524
katowice 9452.72 61.95 2.83 3.05 15.83 0.53 0.43 1979.86 50.26 19.01 337
lodz 8684.23 53.36 2.49 3.94 15.41 0.36 0.27 1975.47 51.76 19.46 1171
bydgoszcz 8333.02 56.05 2.65 2.06 18.51 0.26 0.17 1972.13 53.12 18.01 626
czestochowa 7372.18 55.97 2.56 2.80 10.57 0.21 0.14 1992.08 50.82 19.11 135
radom 6938.81 56.32 2.63 1.90 20.96 0.22 0.19 1977.72 51.40 21.16 150

3. Wizualizacja Zjawiska

3.1 Mapa Statyczna – Średnia Cena za m²

Mapa 15 największych miast Polski – średnia cena mieszkania za m²

Mapa 15 największych miast Polski – średnia cena mieszkania za m²

3.2 Interaktywna Mapa Ofert

3.3 Wykres Słupkowy – Ranking Miast

Wnioski z wizualizacji:

Bezapelacyjna dominacja Warszawy: Stolica wyraźnie dominuje w zestawieniu, osiągając najwyższą średnią cenę na poziomie ok. 18,5 tys. zł/m². Jest również zdecydowanym liderem pod względem wielkości rynku i podaży – posiada ponad dwukrotnie więcej ogłoszeń (6 140) niż drugi w kolejności Kraków (2 712), co doskonale obrazuje wielkość punktu na mapie.

Wyraźny podział na rynki: Analiza ukazuje silną polaryzację. Wyraźnie odznaczają się najdroższe rynki tier-1, do których oprócz Warszawy należą Kraków (17,2 tys. zł/m²) oraz Gdańsk (15,5 tys. zł/m²). Razem z Gdynią i Wrocławiem tworzą one czołówkę miast ze średnią ceną powyżej 13,5 tys. zł/m². Na drugim biegunie znajdują się miasta takie jak Radom (6,9 tys. zł/m²) czy Częstochowa (7,4 tys. zł/m²), gdzie mieszkania są ponad 2,5-krotnie tańsze niż w stolicy.

Geografia i autokorelacja przestrzenna: Najdroższe miasta nie skupiają się w jednej części Polski, lecz stanowią silne, regionalne bieguny wzrostu. Zjawisko wpływu przestrzennego na ceny w miastach ościennych (autokorelacja) można wyraźnie zauważyć jedynie na przykładzie Trójmiasta. W przypadku pozostałych głównych rynków odległości geograficzne między nimi są zbyt duże, by bezpośrednio oddziaływały na siebie nawzajem.


4. Macierz Wag Przestrzennych \(W\)

Uzasadnienie wyboru macierzy odległości odwrotnej:

Analizujemy 15 obserwacji punktowych (jedno miasto = jeden punkt GPS). Klasyczna macierz sąsiedztwa (queen/rook contiguity) jest zdefiniowana dla obszarów z granicami – dla danych punktowych daje wyniki niestabilne.

Macierz oparta na odwrotności odległości euklidesowej (\(w_{ij} = 1/d_{ij}\)) jest metodologicznie poprawna, ponieważ:

  1. Każde miasto ma niezerowe wagi do wszystkich pozostałych – brak arbitralnych progów sąsiedztwa.
  2. Wpływ maleje monotonicznie ze wzrostem odległości (Toblerowe prawo geografii).
  3. Po standaryzacji wierszowej suma wag = 1, co ułatwia interpretację parametrów \(\rho\) (SAR) i \(\lambda\) (SEM).

Zakres odległości między miastami: od 19 km do 638 km.

Macierz wag W (standaryzacja wierszowa) – 9×9
warszawa krakow gdansk gdynia wroclaw poznan rzeszow bialystok lublin
warszawa 0.0000 0.0581 0.0512 0.0482 0.0484 0.0523 0.0576 0.0829 0.0957
krakow 0.0580 0.0000 0.0301 0.0290 0.0616 0.0435 0.0993 0.0359 0.0647
gdansk 0.0370 0.0217 0.0000 0.5689 0.0280 0.0432 0.0198 0.0321 0.0242
gdynia 0.0355 0.0214 0.5815 0.0000 0.0276 0.0420 0.0195 0.0316 0.0238
wroclaw 0.0591 0.0753 0.0474 0.0456 0.0000 0.1235 0.0480 0.0375 0.0462
poznan 0.0604 0.0504 0.0692 0.0658 0.1169 0.0000 0.0382 0.0395 0.0413
rzeszow 0.0754 0.1302 0.0359 0.0347 0.0515 0.0433 0.0000 0.0541 0.1369
bialystok 0.1310 0.0568 0.0702 0.0676 0.0486 0.0540 0.0653 0.0000 0.1077
lublin 0.1107 0.0751 0.0388 0.0373 0.0439 0.0414 0.1211 0.0788 0.0000

5. Badanie Autokorelacji Przestrzennej

5.1 Globalna Statystyka I Morana

Wyniki testu globalnego I Morana
Statystyka Wartość
Moran I statistic I Morana -0.068300
Expectation Oczekiwana wartość E(I) -0.071400
Variance Wariancja Var(I) 0.006908
p-value 0.969800

Interpretacja I Morana:

Statystyka I Morana wynosi -0.0683 (p-value = 0.9698). Wynik ten jest zdecydowanie nieistotny statystycznie. Brak przesłanek do odrzucenia hipotezy zerowej, która zakłada losowy rozkład przestrzenny badanej zmiennej.

Nie mamy potwierdzenia autokorelacji przestrzennej badanego zjawiska, przynajmniej na próbie piętnastu największych miast w Polsce.Oznacza to, że z perspektywy tego testu rozmieszczenie tanich i drogich rynków w kraju jest losowe. Choć możemy obserwować lokalne podobieństwa (np. wspomniane wcześniej wysokie ceny w Trójmieście), to na poziomie całego kraju położenie geograficzne jednego dużego miasta względem drugiego nie determinuje w sposób statystycznie istotny podobieństwa ich cen. Rynki te funkcjonują jako stosunkowo niezależne od siebie “wyspy” cenowe.

5.2 Wykres Rozproszenia Morana (Moran Scatterplot)

Wykres rozproszenia Morana – cena/m² vs. przestrzenny lag ceny

Wykres rozproszenia Morana – cena/m² vs. przestrzenny lag ceny

Interpretacja wykresu Morana:

Interpretacja wykresu rozrzutu Morana: Każdy punkt na wykresie reprezentuje jedno miasto. Oś X pokazuje standaryzowaną cenę mieszkania, natomiast oś Y przedstawia przestrzenną średnią ważoną cen sąsiadów (tzw. wygładzenie przestrzenne). Analiza pozwala zidentyfikować dwa klastry oraz dwa typy punktów nietypowych (outlierów):

  • Klaster HH (wysoka–wysoka / czerwony): Do tego klastra należą Gdańsk oraz Gdynia. Potwierdza to silną, lokalną relację w obrębie Trójmiasta – są to drogie rynki położone w swoim bezpośrednim sąsiedztwie.
  • Klaster LL (niska–niska / niebieski): W tym kwadrancie znajdują się Lublin oraz Rzeszów, reprezentujące rynki o niższych cenach, które są otoczone przez regiony o podobnie niższym profilu cenowym.
  • Outlierzy HL (wysoka–niska / pomarańczowy): Wyraźnymi nietypowymi punktami na mapie są Warszawa, Kraków, Wrocław oraz Poznań. Są to dynamiczne, drogie rynki, które jednak w skali geograficznej sąsiadują z rynkami o znacznie niższych średnich cenach.
  • Outlierzy LH (niska–wysoka / zielony): Do tej grupy należy najwięcej miast: Radom, Bydgoszcz, Częstochowa, Łódź, Katowice, Białystok oraz Szczecin. Są to ośrodki o relatywnie niskich cenach własnych, zlokalizowane w sąsiedztwie droższych rynków.

Główny wniosek: Dominująca liczba miast w Polsce to outlierzy przestrzeni (HL oraz LH). Uniemożliwia to sformułowanie jednej, ogólnokrajowej reguły dotyczącej klastrowania cen mieszkań. Zjawisko to jest bezpośrednim efektem geograficznego rozproszenia największych, najdroższych metropolii po całej mapie kraju.

5.3 Lokalna Statystyka I Morana (LISA)

Mapa LISA – lokalna autokorelacja przestrzenna cen mieszkań

Mapa LISA – lokalna autokorelacja przestrzenna cen mieszkań

Lokalne statystyki I Morana (LISA) dla 15 miast
Miasto Śr. cena/m² Li (Ii) p-value Klaster Istotny (α=0.05)
warszawa 18518 -0.7388 0.0523 Nieistotny ❌ Nie |
krakow 17155 -0.5945 0.1713 Nieistotny ❌ Nie |
gdynia 14108 0.4495 0.2515 Nieistotny ❌ Nie |
gdansk 15451 0.3981 0.4211 Nieistotny ❌ Nie |
radom 6939 -0.1990 0.6870 Nieistotny ❌ Nie |
bydgoszcz 8333 -0.1428 0.4554 Nieistotny ❌ Nie |
wroclaw 13649 -0.1047 0.2767 Nieistotny ❌ Nie |
bialystok 10213 -0.0345 0.4898 Nieistotny ❌ Nie |
szczecin 9703 -0.0317 0.8272 Nieistotny ❌ Nie |
czestochowa 7372 -0.0289 0.7237 Nieistotny ❌ Nie |
lublin 10065 0.0248 0.5607 Nieistotny ❌ Nie |
katowice 9453 -0.0154 0.9409 Nieistotny ❌ Nie |
lodz 8684 -0.0068 0.6336 Nieistotny ❌ Nie |
rzeszow 10935 0.0047 0.7693 Nieistotny ❌ Nie |
poznan 11661 -0.0043 0.5367 Nieistotny ❌ Nie |

Interpretacja mapy klastrów LISA:

Test LISA potwierdził wnioski z globalnego testu Morana – analizując 15 największych miast w Polsce, nie mamy do czynienia z istotną lokalną korelacją przestrzenną. Najbliżej istotności statystycznej na poziomie 0,05 była Warszawa (p-value = 0,0523). Gdyby przyjąć poziom istotności 0,1, Warszawa stałaby się outlierem przestrzennym (typu HL) w naszym modelu – miastem o wysokich cenach, podczas gdy średnia cena jego sąsiadów jest znacznie niższa. Przy przyjętym zbiorze danych jest to w pełni logiczne, gdyż jednym z najbliższych sąsiadów Warszawy w sieci powiązań jest Radom, stanowiący najtańsze z 15 badanych miast.


6. Modelowanie Ekonometryczne

6.1 Przygotowanie Zmiennych

Specyfikacja modelu bazowego:

\[\ln(\text{cena/m}^2_i) = \beta_0 + \beta_1 \ln(\text{metraż}_i) + \beta_2 \text{pokoje}_i + \beta_3 \ln(\text{odl. centrum}_i) + \beta_4 \text{POI}_i + \beta_5 \text{winda}_i + \beta_6 \text{parking}_i + \beta_7 \text{rok budowy}_i + \varepsilon_i\]

Zmienne transformowane logarytmicznie: avg_price_m2, avg_sqm, avg_centre_dist – stabilizacja wariancji i interpretacja elastyczności procentowych.

Liczba miast w modelu: 15

6.2 Model OLS – Punkt Odniesienia

## 
## Call:
## lm(formula = formula_ols, data = df_miasta)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.09888 -0.07271 -0.01516  0.05656  0.15974 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)  
## (Intercept)      -11.744192  15.144910  -0.775   0.4635  
## log_sqm           -0.279108   1.908614  -0.146   0.8879  
## avg_rooms          0.241949   0.820026   0.295   0.7765  
## log_centre_dist    0.418104   0.239462   1.746   0.1243  
## avg_poi_count      0.027841   0.012511   2.225   0.0614 .
## avg_has_elevator   0.788483   0.646860   1.219   0.2623  
## avg_has_parking   -0.299248   0.460631  -0.650   0.5366  
## avg_build_year     0.010217   0.006425   1.590   0.1558  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1289 on 7 degrees of freedom
## Multiple R-squared:  0.9072, Adjusted R-squared:  0.8144 
## F-statistic: 9.777 on 7 and 7 DF,  p-value: 0.00376
Wykresy diagnostyczne modelu OLS

Wykresy diagnostyczne modelu OLS

Interpretacja modelu OLS:

Istotność zmiennych

Na poziomie istotności 0,05 żadna zmienna modelu nie jest statystycznie istotna. Głównym powodem powyższego może być bardzo mała próba (15 miast) przy stosunkowo wielu zmiennych, co wpływa na stopnie swobody.

Dopasowanie modelu

Współczynnik R2 wynosi aż 0,9072, więc zawarte w modelu zmienne wyjaśniają aż 90,7% zmienności cen mieszkań. Dodatkowo należy zaznaczyć, że pomimo braku istotności zmiennych pojedynczo, model jako całość (test F) jest istotny (p-value 0,00376).

Zgodność współczynników z teorią

Kierunki zmian w większości są logiczne i zgodne z teorią ekonomiczną. Wyjątki stanowią odległość od centrum z dodatnim znakiem oraz posiadanie parkingu ze znakiem ujemnym. Wynika to najprawdopodobniej ze współliniowości danych. Np. dla parkingu ten brak logiki można objasnić problemem z miejscami parkingowymi który ma miejsce w dużych aglomeracjach (pamiętamy, że model nie analizował każdego mieszkania osobno, a wyciągał średnią wartość dla każdego miasta).

Diagnostyka reszt modelu

Normalność reszt: Na wykresie Q-Q punkty trzymają się blisko linii prostej, więc założenie o normalności rozkładu jest w miarę ok (wyskakują tylko pojedyncze obserwacje: 2,4 i 13)

Homoskedastyczność: Wykresy reszt lekko falują, co pokazuje, że wariancja nie jest idealnie równa (częsty problem przy porównywaniu skrajnie różnych wielkościowo miast).

Wartości wpływowe: Wykres odległości Cooka (Cook’s distance) wyraźnie pokazuje, że miasta o numerach 12, 14 i 4 mocno odstają i mocno “ciągną” za sobą wyniki całego modelu. Gdybyśmy je wyrzucili, wyniki mogłyby wyglądać zupełnie inaczej.

6.3 Testy ex-ante – Mnożniki Lagrange’a (LM)

Testy Mnożników Lagrange’a (diagnostyka ex-ante)
Test Statystyka p_value Wynik
LMlag (SAR) 0.3230 0.5698 ❌ Nieistotny |
LMerr (SEM) 1.5147 0.2184 ❌ Nieistotny |
RLMlag (robust SAR) 5.3919 0.0202 ✅ Istotny |
RLMerr (robust SEM) 6.5836 0.0103 ✅ Istotny |
SARMA 6.9066 0.0316 ✅ Istotny |
Schemat wyboru modelu na podstawie testów LM:

Analiza wyników testów

Przeprowadzone testy mnożników Lagrange’a (LM) dają niejednoznaczny obraz przy klasycznym podejściu, ponieważ podstawowe testy LMlag (p-value = 0,5698) oraz LMerr (p-value = 0,2184) okazały się statystycznie nieistotne.

Jednak w przypadku małej próby (15 miast) kluczowe są wersje odporne. Zarówno RLMlag (robust SAR) z p-value = 0,0202, jak i RLMerr (robust SEM) z p-value = 0,0103 są istotne statystycznie na poziomie alpha = 0,05. Dodatkowo test SARMA potwierdza obecność efektów przestrzennych (p-value = 0,0316).

Ponieważ w wersjach odpornych (robust) oba testy okazały się istotne, zgodnie z procedurą Anselina konieczne jest porównanie wartości statystyk testowych dla RLMlag oraz RLMerr:

Statystyka RLMlag: 5.3919

Statystyka RLMerr: 6.5836

Z racji tego, że wartość statystyki dla błędu przestrzennego jest większa (RLMerr > RLMlag), ostatecznym i najbardziej wskazanym modelem do dalszej analizy jest model błędu przestrzennego – SEM (Spatial Error Model).

Wybór ten oznacza, że przestrzenna zależność na rynku mieszkaniowym w Polsce nie wynika bezpośrednio z faktu, że ceny w jednym mieście „zarażają” sąsiadów, ale ukrywa się w czynnikach losowych (składniku resztowym) – takich jak np. regionalne uwarunkowania ekonomiczne, specyfika lokalnego popytu czy położenie geograficzne, których nie ujęto bezpośrednio w zmiennych modelu bazowego OLS.

6.4 Model SAR (Spatial Autoregressive / Spatial Lag)

\[Y = \rho W Y + X\beta + \varepsilon, \quad \varepsilon \sim \mathcal{N}(0, \sigma^2 I)\]

## 
## Call:lagsarlm(formula = formula_ols, data = df_miasta, listw = W_listw, 
##     method = "eigen")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.099374 -0.063886 -0.024942  0.043445  0.166292 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##                     Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)      -15.7147886  11.1327478 -1.4116 0.1580730
## log_sqm           -0.5314883   1.2825461 -0.4144 0.6785805
## avg_rooms          0.2808945   0.5474201  0.5131 0.6078645
## log_centre_dist    0.4442296   0.1610898  2.7577 0.0058218
## avg_poi_count      0.0296359   0.0085716  3.4574 0.0005454
## avg_has_elevator   0.7063371   0.4356245  1.6214 0.1049242
## avg_has_parking   -0.1168087   0.3518105 -0.3320 0.7398729
## avg_build_year     0.0111451   0.0044052  2.5300 0.0114063
## 
## Rho: 0.31895, LR test value: 0.48142, p-value: 0.48778
## Asymptotic standard error: 0.34386
##     z-value: 0.92757, p-value: 0.35363
## Wald statistic: 0.86039, p-value: 0.35363
## 
## Log likelihood: 15.4025 for lag model
## ML residual variance (sigma squared): 0.0074041, (sigma: 0.086047)
## Nagelkerke pseudo-R-squared: 0.91014 
## Number of observations: 15 
## Number of parameters estimated: 10 
## AIC: -10.805, (AIC for lm: -12.324)
## LM test for residual autocorrelation
## test value: 6.8516, p-value: 0.0088564

6.5 Model SEM (Spatial Error Model)

\[Y = X\beta + u, \quad u = \lambda W u + \varepsilon, \quad \varepsilon \sim \mathcal{N}(0, \sigma^2 I)\]

## 
## Call:errorsarlm(formula = formula_ols, data = df_miasta, listw = W_listw, 
##     method = "eigen")
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.1162423 -0.0497783 -0.0061662  0.0569668  0.1316945 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                     Estimate  Std. Error z value Pr(>|z|)
## (Intercept)      -13.9036072  10.4556559 -1.3298 0.183594
## log_sqm            0.5082878   1.2536328  0.4055 0.685145
## avg_rooms         -0.0691244   0.5394511 -0.1281 0.898039
## log_centre_dist    0.3623553   0.1666084  2.1749 0.029638
## avg_poi_count      0.0274783   0.0081207  3.3837 0.000715
## avg_has_elevator   1.0011919   0.4345297  2.3041 0.021218
## avg_has_parking   -0.4426881   0.2534080 -1.7469 0.080648
## avg_build_year     0.0101269   0.0044230  2.2896 0.022046
## 
## Lambda: -1.054, LR test value: 2.8209, p-value: 0.093047
## Asymptotic standard error: 0.44073
##     z-value: -2.3914, p-value: 0.016784
## Wald statistic: 5.7188, p-value: 0.016784
## 
## Log likelihood: 16.57221 for error model
## ML residual variance (sigma squared): 0.0056932, (sigma: 0.075453)
## Nagelkerke pseudo-R-squared: 0.92312 
## Number of observations: 15 
## Number of parameters estimated: 10 
## AIC: -13.144, (AIC for lm: -12.324)

6.6 Model SDM (Spatial Durbin Model)

\[Y = \rho W Y + X\beta + WX\theta + \varepsilon\]

## 
## Call:lagsarlm(formula = formula_sdm, data = df_miasta, listw = W_listw, 
##     type = "Durbin", method = "eigen")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.284389 -0.089857 -0.020955  0.126230  0.266344 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                      Estimate Std. Error z value Pr(>|z|)
## (Intercept)         -9.649547  31.983487 -0.3017 0.762878
## log_sqm              0.475415   1.066335  0.4458 0.655712
## log_centre_dist      0.596719   0.185176  3.2224 0.001271
## avg_poi_count        0.016840   0.015351  1.0970 0.272630
## lag.log_sqm          5.033304   6.945206  0.7247 0.468626
## lag.log_centre_dist  0.081195   1.248625  0.0650 0.948152
## lag.avg_poi_count    0.032391   0.113439  0.2855 0.775235
## 
## Rho: -0.54296, LR test value: 0.46661, p-value: 0.49455
## Asymptotic standard error: 0.57275
##     z-value: -0.94798, p-value: 0.34314
## Wald statistic: 0.89867, p-value: 0.34314
## 
## Log likelihood: 5.655266 for mixed model
## ML residual variance (sigma squared): 0.026681, (sigma: 0.16334)
## Nagelkerke pseudo-R-squared: 0.67041 
## Number of observations: 15 
## Number of parameters estimated: 9 
## AIC: 6.6895, (AIC for lm: 5.1561)
## LM test for residual autocorrelation
## test value: 3.4998, p-value: 0.061377

7. Wybór Najlepszego Modelu

7.1 Porównanie AIC

Porównanie modeli wg kryterium AIC (niższe = lepsze)
Model LogLik AIC Parametry Delta_AIC Preferowany
OLS OLS 15.162 -12.324 9 0.820
SAR SAR 15.402 -10.805 10 2.339
SEM SEM 16.572 -13.144 10 0.000 ⭐ Najlepszy |
SDM SDM 5.655 6.689 9 19.833

7.2 Testy LR (Likelihood Ratio)

## ── LR test: OLS vs SAR ────────────────────────────────
## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = 0.48142, df = 1, p-value = 0.4878
## sample estimates:
## Log likelihood of model_sar Log likelihood of model_ols 
##                    15.40250                    15.16179
## 
## ── LR test: OLS vs SEM ────────────────────────────────
## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = 2.8209, df = 1, p-value = 0.09305
## sample estimates:
## Log likelihood of model_sem Log likelihood of model_ols 
##                    16.57221                    15.16179
## 
## ── LR test: SAR vs SDM ────────────────────────────────
## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = -19.494, df = 1, p-value = 1.009e-05
## sample estimates:
## Log likelihood of model_sdm Log likelihood of model_sar 
##                    5.655266                   15.402496 
## 
## 
## ── LR test: SEM vs SDM ────────────────────────────────
## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = -21.834, df = 1, p-value = 2.973e-06
## sample estimates:
## Log likelihood of model_sdm Log likelihood of model_sem 
##                    5.655266                   16.572215

7.3 Schemat Wyboru Modelu

Parametry przestrzenne oszacowanych modeli
Model Parametr Wartość Błąd std. Z-stat p-value Istotność
rho SAR ρ (rho) 0.319 0.3439 0.9276 0.3536 ❌ Nieistotny |
lambda SEM λ (lambda) -1.054 0.4407 -2.3916 0.0168 ✅ p < 0.05 |
rho1 SDM ρ (rho) -0.543 0.5727 -0.9481 0.3431 ❌ Nieistotny |
Wybór najlepszego modelu – uzasadnienie:

Tak jak przedstawiono w punkcie 6.3, testy Ex-Ante Lagrange’a (LM) wskazały na zasadność zastosowania modelu SEM (RLM err > RLM lag).

Oszacowanie parametru strukturalnego jednoznacznie potwierdza ten kierunek. Dla modeli SAR oraz SDM parametry ρ okazały się statystycznie nieistotne. Jedynym modelem, w którym wpływ przestrzenny jest istotny statystycznie, jest model SEM.

Parametr λ dla modelu SEM wynosi -1.054, a jego p-value = 0.0168, co oznacza pełną istotność statystyczną na poziomie α = 0.05.

Wybranym modelem końcowym jest więc model SEM.

Ujemna wartość parametru λ przy istotnym statystycznie błędzie przestrzennym oznacza, że niewyjaśnione reszty modelu mają tendencję do lokalnego różnicowania się w przestrzeni.

Oznacza to, że — pokazując na przykładzie Warszawy i Radomia — najdroższa Warszawa jest najbliższym sąsiadem dla jednego z najtańszych rynków mieszkaniowych. Model interpretuje to jako sytuację, w której największe i najdroższe miasta „zasysają” zasoby z pozostałych okolicznych dużych ośrodków miejskich.

Dzięki wykryciu tej zależności model SEM zapewnia bardziej stabilne i wiarygodne oszacowania parametrów niż klasyczny model KMNK (OLS).

7.4 Tabela Zbiorcza – stargazer

Porównanie wyników modeli – OLS (spatial modele poniżej)
Dependent variable:
log(Śr. cena/m²)
log(Metraż) -0.279
(1.909)
Śr. l. pokoi 0.242
(0.820)
log(Odl. centrum) 0.418
(0.239)
Śr. POI 0.028*
(0.013)
Śr. winda 0.788
(0.647)
Śr. parking -0.299
(0.461)
Śr. rok budowy 0.010
(0.006)
Constant -11.744
(15.145)
Observations 15
R2 0.907
Adjusted R2 0.814
Note: p<0.1; p<0.05; p<0.01
Model SAR: ρ = 0.319 | Model SEM: λ = -1.054 | Model SDM: ρ = -0.543

Zwykły model OLS poległ: Mimo że ogólny wynik R2 jest wysoki (0,907), to przez małą próbę (15 miast) model OLS nie wykazał istotności statystycznej dla prawie żadnej zmiennej. Działał tak, jakby te miasta były od siebie całkowicie odizolowane.

Geografia ma znaczenie: W rzeczywistości miasta na siebie wpływają, co potwierdza parametr λ = -1,054 z modelu SEM, który jako jedyny wyszedł w pełni istotny (p < 0,05).

Minus oznacza konkurencję: Ujemna lambda pokazuje, że polskie miasta nie współpracują, tylko konkurują ze sobą o kapitał, studentów i inwestorów. Jeśli jakieś ukryte zjawisko (np. prestiż stolicy) sztucznie podbija ceny w Warszawie, to odbywa się to kosztem wysysania zasobów z sąsiednich miast (np. Radomia), przez co tam ceny stają się relatywnie niższe.


8. Efekty Bezpośrednie, Pośrednie i Całkowite

## === Efekty bezpośrednie, pośrednie i całkowite – Model SAR ===
## Impact measures (lag, exact):
##                             Direct    Indirect       Total
## log_sqm dy/dx          -0.53969977 -0.24069851 -0.78039828
## avg_rooms dy/dx         0.28523426  0.12721047  0.41244474
## log_centre_dist dy/dx   0.45109287  0.20118108  0.65227396
## avg_poi_count dy/dx     0.03009374  0.01342138  0.04351512
## avg_has_elevator dy/dx  0.71724995  0.31988340  1.03713335
## avg_has_parking dy/dx  -0.11861338 -0.05289990 -0.17151328
## avg_build_year dy/dx    0.01131733  0.00504737  0.01636470
## ========================================================
## Simulation results ( variance matrix):
## ========================================================
## Simulated standard errors
##                             Direct   Indirect      Total
## log_sqm dy/dx          1.418692589 4.72577768 5.58872684
## avg_rooms dy/dx        0.605351416 2.07716094 2.43805915
## log_centre_dist dy/dx  0.200077799 1.54058835 1.67335233
## avg_poi_count dy/dx    0.014646310 0.14805545 0.16009324
## avg_has_elevator dy/dx 0.610889655 5.11247278 5.55053612
## avg_has_parking dy/dx  0.374376307 1.16864466 1.39815376
## avg_build_year dy/dx   0.006503999 0.05852952 0.06351532
## 
## Simulated z-values:
##                            Direct    Indirect       Total
## log_sqm dy/dx          -0.3479612 -0.14729932 -0.21288459
## avg_rooms dy/dx         0.4349621  0.10333296  0.19603467
## log_centre_dist dy/dx   2.3277145  0.34816949  0.59886363
## avg_poi_count dy/dx     2.1671090  0.27931534  0.45657335
## avg_has_elevator dy/dx  1.3267448  0.20105653  0.33120958
## avg_has_parking dy/dx  -0.4163969  0.05735999 -0.06355215
## avg_build_year dy/dx    1.8513429  0.27577567  0.44370633
## 
## Simulated p-values:
##                        Direct   Indirect Total  
## log_sqm dy/dx          0.727869 0.88290  0.83142
## avg_rooms dy/dx        0.663590 0.91770  0.84458
## log_centre_dist dy/dx  0.019927 0.72771  0.54926
## avg_poi_count dy/dx    0.030227 0.78000  0.64798
## avg_has_elevator dy/dx 0.184593 0.84065  0.74049
## avg_has_parking dy/dx  0.677120 0.95426  0.94933
## avg_build_year dy/dx   0.064120 0.78272  0.65725
Efekty bezpośrednie, pośrednie i całkowite – Model SAR
Zmienna Efekt bezpośredni Efekt pośredni (spillover) Efekt całkowity
log_sqm dy/dx -0.5397 -0.2407 -0.7804
avg_rooms dy/dx 0.2852 0.1272 0.4124
log_centre_dist dy/dx 0.4511 0.2012 0.6523
avg_poi_count dy/dx 0.0301 0.0134 0.0435
avg_has_elevator dy/dx 0.7172 0.3199 1.0371
avg_has_parking dy/dx -0.1186 -0.0529 -0.1715
avg_build_year dy/dx 0.0113 0.0050 0.0164

Definicja efektów (Model SAR):

Co oznaczają poszczególne efekty:

Efekt bezpośredni: Jak zmiana cechy mieszkania (np. metrażu) w danym mieście wpływa na cenę w tym samym mieście.

Efekt pośredni (spillover): Efekt domina. Jak zmiana w jednym mieście „rozlewa się” i zmienia ceny u jego sąsiadów.

Efekt całkowity: Łączna siła uderzenia (suma efektu w mieście i u sąsiadów).

Kluczowe wnioski z tabeli:

Największy wpływ na ceny ma winda. Jej efekt całkowity wynosi aż 1.0371. Jeśli w mieście przybywa bloków z windami, ceny skaczą mocno u nas (0.7172), ale impuls podbija też ceny w miastach sąsiednich (0.3199).

Metraż zgodnie z logiką ma minus (-0.7804). Im większe mieszkania buduje się na rynku, tym niższa cena za metr u nas (-0.5397) oraz u sąsiadów (-0.2407).

Dziwne zachowanie odległości od centrum i parkingów. Podobnie jak w OLS, te zmienne mają odwrócone znaki (parking ma minus, centrum ma plus). To anomalie wynikające ze specyfiki próby 15 miast, ale model SAR pokazuje, że te błędy również mocno “zarażają” sąsiednie rynki.

Interpretacja efektów: Wzrost średniego metrażu o 1% w danym mieście powoduje spadek ceny za m^2 o 0.54% w tym samym mieście (efekt bezpośredni) oraz dodatkowy spadek o 0.24% w miastach sąsiednich (efekt pośredni). Łączny ujemny wpływ tej zmiennej na rynek wynosi -0.78%. Z kolei zmienną o najsilniejszym, dodatnim wpływie całkowitym jest obecność windy (efekt łączny: 1.0371), która generuje silny impuls wzrostowy zarówno lokalnie, jak i na rynkach ościennych.


9. Diagnostyka Reszt i Wizualizacja Wyników

9.1 Mapa Reszt Modelu OLS

Mapa reszt modelu OLS – identyfikacja outlierów przestrzennych

Mapa reszt modelu OLS – identyfikacja outlierów przestrzennych

Geograficzny rozkład błędów modelu OLS wyraźnie potwierdza występowanie autokorelacji przestrzennej składnika losowego. Model OLS systematycznie i silnie niedoszacowuje cen nieruchomości w kluczowych rynkach regionalnych (czerwone punkty: Kraków ze skrajną resztą 0,16, Gdynia 0,152 oraz Bydgoszcz 0,132. Jednocześnie model przeszacowuje rynki leżące w strefach oddziaływania liderów cenowych (niebieskie punkty: Warszawa -0,078, Radom -0,068, Lublin -0,099 czy Katowice -0,055. Taki układ przestrzenny (sąsiadowanie skrajnie dodatnich i ujemnych błędów) ostatecznie uzasadnia konieczność zastosowania modelu błędu przestrzennego (SEM).

9.2 Wartości Dopasowane vs Rzeczywiste

Wartości dopasowane vs rzeczywiste – OLS, SAR, SEM

Wartości dopasowane vs rzeczywiste – OLS, SAR, SEM

Wykresy porównawcze idealnie obrazują przewagę modeli przestrzennych nad klasycznym OLS. W modelu OLS kropki reprezentujące miasta mocno odlatują od przerywanej linii idealnego dopasowania (y=x), generując duże dysproporcje (szczególnie dla Krakowa, Warszawy i Gdyni). Model SAR koryguje te błędy tylko nieznacznie. Dopiero model SEM dokonuje matematycznej korekty – dzięki uwzględnieniu parametru λ kropki niemal idealnie układają się wzdłuż linii skośnej. Oznacza to, że model SEM skutecznie wychwycił lokalną rywalizację miast o kapitał, eliminując systematyczny błąd prognozy i dając najbardziej zbliżone do rzeczywistości, wiarygodne wyceny.


10. Podsumowanie i Wnioski

10.1 Tabela Zbiorcza Wyników

Zbiorcze podsumowanie wyników analizy przestrzennej
Kryterium Wartość
I Morana (cena/m²) -0.0683
p-value I Morana 0.9698
LMlag (SAR) 0.3230
p-value LMlag 0.5698
LMerr (SEM) 1.5147
p-value LMerr 0.2184
AIC – OLS -12.3200
AIC – SAR -10.8000
AIC – SEM -13.1400
AIC – SDM 6.6900
ρ (SAR) 0.3190
λ (SEM) -1.0540

10.2 Wnioski Ekonomiczne

1. Autokorelacja przestrzenna cen mieszkań

Statystyka I Morana wynosi -0.0683 i jest nieistotna statystycznie (p-value = 0.9698). Oznacza to, że ceny mieszkań w polskich miastach są rozłożone losowo w przestrzeni, jeśli patrzymy na nie bezpośrednio, bez uwzględnienia innych zmiennych. Analiza LISA potwierdziła ten wynik – dla żadnego z 15 miast nie udało się wyodrębnić statystycznie istotnego klastra typu HH czy LL na poziomie alpha = 0.05.

2. Wybór modelu przestrzennego

Na podstawie testów LM i kryterium AIC jako optymalny model wybrany został SEM (Spatial Error Model). Choć tradycyjne testy LM były nieistotne, to odporny test RLMerr (robust SEM) okazał sięistotny (p = 0.0103) i przewyższył statystyką RLMlag (6.58 > 5.39). Model SEM charakteryzuje się również najniższą wartością kryterium informacyjnego (AIC = -13.14), co potwierdza jego najwyższą jakość.

3. Efekty bezpośrednie i pośrednie

Wzrost średniego metrażu o 1% w danym mieście powoduje spadek ceny za m^2 o 0.54% w tym samym mieście (efekt bezpośredni) oraz dodatkowy spadek (efekt spillover) o 0.24% w pozostałych miastach (efekt pośredni). Łączny wpływ całkowity wynosi -0.78%.

4. Parametr ρ (SAR) / λ (SEM)

Parametr lambda wynosi -1.0540. Jego ujemny kierunek wskazuje na przestrzenną rywalizację i kompensację błędów. Losowe czynniki, które stymulują wzrost cen w jednym silnym ośrodku (np. prestiż Warszawy czy Krakowa), działają jak magnes i wysysają kapitał inwestycyjny z sąsiednich rynków, relatywnie zaniżając ich ceny.

5. Główne determinanty cen

Spośród zmiennych objaśniających największy wpływ na cenę za m^2 mają: winda (beta = 0.788) oraz liczba punktów POI (beta = 0.028, jedyna zmienna wykazująca cechy istotności w OLS). Logika ekonomiczna wskazuje, że nowoczesna infrastruktura budynku (winda) drastycznie podnosi standard i wycenę lokalu, a bliskość punktów użyteczności publicznej (POI) generuje premię za lokalizację, kluczową dla rynków w największych metropoliach.

10.3 Ograniczenia Badania

  1. Mała próba (N = 15): Analiza na poziomie miast opiera się na zaledwie 15 obserwacjach – ograniczona moc statystyczna testów.

  2. Agregacja danych: Uśrednienie cen do poziomu miasta pomija zróżnicowanie wewnątrzmiejskie (centrum vs peryferia).

  3. Brak zmiennych makroekonomicznych: Model nie uwzględnia stóp procentowych, dochodów per capita ani bezrobocia.

  4. Statyczność modelu: Dane z jednego miesiąca (VI 2024) – brak analizy zmian w czasie.

  5. Specyfika geograficzna i dystans między metropoliami: Nasza próba obejmuje tylko 15 największych miast w skali całego kraju. W polskiej rzeczywistości geograficznej te ośrodki (poza nielicznymi wyjątkami jak Trójmiasto) są od siebie oddalone o setki kilometrów. Przy tak dużych odległościach fizycznych bezpośrednie interakcje przestrzenne na rynku nieruchomości są naturalnie słabsze i trudniejsze do uchwycenia przez klasyczne macierze sąsiedztwa, co mogło wpłynąć na brak istotności testów globalnych (np. I Morana).


Kurs: Ekonometria Przestrzenna – AG II N, semestr II | Prowadzący: dr Aleksandra Kordalska