Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update time_conversion.py #3

Merged

Conversation

steveshaoucsb
Copy link
Contributor

Add adjustments in handle_endpoint_differences to resolve issues related to quarterly data:
Setting C_mask as "None" and immediately return it so that "C_mask = C_mask or np.full(nl, True)" in ts_disagg.py/build_conversion_matrix will do "np.full(nl, True)" to make the matrix size match up.

Note: Not sure it is absolutely correct but it works!

Add adjustments in handle_endpoint_differences to resolve issues related to quarterly data
@jessegrabowski
Copy link
Owner

jessegrabowski commented Aug 6, 2024

This looks really great! All the old tests are passing, so it should be good to merge. Can you just add a regression test? Basically add test_chow_lin_Q_to_M to the DisaggregationTests in test_disaggregation.py, and run exactly what was failing for you to make sure its fixed now. You can convert imports_q to monthly using exports_m as an indicator.

Add test files
@jessegrabowski
Copy link
Owner

No need, you can keep working in this PR.

You can check this test for what I'm currently doing. The file you uploaded is the "ground truth" from the R function?

@steveshaoucsb
Copy link
Contributor Author

steveshaoucsb commented Aug 6, 2024

Yes, I did run it in R on my side. Here is the code for your reference:

library(tidyverse)
library(tsdisagg2)
library(tempdisagg)

data(tempdisagg)
exports_m <- read.csv("exports_m.csv")
imports_q <- read.csv("imports_q.csv")

view(exports_m)
view(imports_q)

# no indicator (Denton-Cholette)
result <- td(imports.q ~ exports.m + 1, to="monthly", method="chow-lin-maxlog")
monthly_disaggregated <- predict(result)


monthly_disaggregated_df <- data.frame(Value = as.numeric(monthly_disaggregated))

view(monthly_disaggregated_df)
write.csv(monthly_disaggregated_df, file = "R_Output_chow-lin_QtoM.csv", row.names = FALSE)

Is this match what you do to generate "ground truth" in R as well?

@steveshaoucsb
Copy link
Contributor Author

steveshaoucsb commented Aug 7, 2024

Ok so I double check the results on my side from both Python and R, and found two completely different chow-lin regression results. Not sure where is the problem. But I post it here so that you, and my teammate who is also working on this, can take a look:

Results from R:

Call:

td(formula = imports.q ~ exports.m + 1, to = "monthly", method = "chow-lin-maxlog")

Residuals:
    Min      1Q  Median      3Q     Max 
-1388.9  -215.1   -60.0   111.7  1489.9 

Coefficients:
           Estimate   Std. Error  t value    Pr(>|t|)    
(Intercept) 88.08168   31.30688   2.813  0.00553 ** 
exports.m    0.52749    0.01055  49.988  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

'chow-lin-maxlog' disaggregation with 'sum' conversion
158 low-freq. obs. converted to 474 high-freq. obs.
Adjusted R-squared: 0.9409	AR1-Parameter: 0.701

Python code used to generate the result:

imports_q = pd.read_csv('imports_q.csv')
exports_m = pd.read_csv('exports_m.csv')

exports_m.index = pd.date_range(start="1972-01-01", freq="MS", periods=exports_m.shape[0])
exports_m.columns = ["exports"]

imports_q.index = pd.date_range(start="1972-01-01", freq="QS-OCT", periods=imports_q.shape[0])
imports_q.columns = ["imports"]

disaggregated_lf = disaggregate_series(
    imports_q,
    exports_m.assign(intercept=1),
    method="chow-lin",
    agg_func="sum",
    optimizer_kwargs={"method": "powell"}
    )

Results from Python:

Dependent Variable: quarterly_imports
GLS Estimates using Chow-Lin's covariance matrix
N = 474		df = 470
Adj r2 = 0.3580

Variable             coef         sd err              t        P > |t|         [0.025         0.975]
----------------------------------------------------------------------------------------------------
exports            0.1093         0.0098        11.1505         0.0000         0.0900         0.1286
intercept       1219.6370       667.4767         1.8272         0.0341       -91.9709      2531.2450

rho                0.9955
sigma.sq        4077.4387

@jessegrabowski
Copy link
Owner

The scale difference in the intercept suggests that either the optimizer isn't converging, or the data isn't the same. I'll have some time to take a look later today and I'll see if I can help get to the bottom of this.

@jessegrabowski
Copy link
Owner

I foolishly didn't add a flag to return the raw optimizer output for debugging, I would suggest doing that as well and checking if it's converging. You can also try using {'method':'Nelder-Mead'}, that should be really robust on low dimensional problems like this one, just to rule out something like that.

@steveshaoucsb
Copy link
Contributor Author

steveshaoucsb commented Aug 7, 2024

Double checked on my side by comparing two dataframe together. The data used in the python side is the same as the R side! {'method':'Nelder-Mead'} didn't yield good results:

results after optimization:        
       message: Optimization terminated successfully.
       success: True
        status: 0
           fun: 1171.0035961399985
             x: [ 1.000e+00  4.149e+03]
           nit: 73
          nfev: 146
 final_simplex: (array([[ 1.000e+00,  4.149e+03],
                       [ 1.000e+00,  4.149e+03],
                       [ 1.000e+00,  4.149e+03]]), array([ 1.171e+03,  1.171e+03,  1.171e+03]))
                       

Dependent Variable: quarterly_imports
GLS Estimates using Chow-Lin's covariance matrix
N = 474		df = 470
Adj r2 = 0.3250

Variable             coef         sd err              t        P > |t|         [0.025         0.975]
----------------------------------------------------------------------------------------------------
exports            0.0980         0.0099         9.9202         0.0000         0.0786         0.1175
intercept       1423.3628     20345.0355         0.0700         0.4721    -38555.1237     41401.8494

rho                1.0000
sigma.sq        4148.9850

@jessegrabowski
Copy link
Owner

Gross. Can you confirm that the feature matrix we build and ultimate use on the RHS of the regression matches R?

@steveshaoucsb
Copy link
Contributor Author

Just checked. Both matrix matches with each other!
For your reference, I put both of them together for clear comparison:

     Intercept  exports.m.Python  Intercept  exports.m.R
0            1        451.701000          1   451.701000
1            1        439.238000          1   439.238000
2            1        541.700000          1   541.700000
3            1        475.392000          1   475.392000
4            1        465.519000          1   465.519000
5            1        515.980000          1   515.980000
6            1        488.318000          1   488.318000
7            1        428.756000          1   428.756000
8            1        425.488000          1   425.488000
9            1        488.820000          1   488.820000
10           1        490.247000          1   490.247000
11           1        560.327000          1   560.327000
12           1        504.986000          1   504.986000
13           1        515.158000          1   515.158000
14           1        515.610000          1   515.610000
15           1        484.418000          1   484.418000
16           1        591.752000          1   591.752000
17           1        502.288000          1   502.288000
18           1        542.665000          1   542.665000
19           1        511.883000          1   511.883000
20           1        520.176000          1   520.176000
21           1        579.106000          1   579.106000
22           1        580.689000          1   580.689000
23           1        492.376000          1   492.376000
24           1        670.063000          1   670.063000
25           1        684.526000          1   684.526000
26           1        693.245000          1   693.245000
27           1        725.774000          1   725.774000
28           1        749.077000          1   749.077000
29           1        643.120000          1   643.120000
30           1        741.455000          1   741.455000
31           1        582.606000          1   582.606000
32           1        601.865000          1   601.865000
33           1        658.485000          1   658.485000
34           1        624.968000          1   624.968000
35           1        514.737000          1   514.737000
36           1        635.638000          1   635.638000
37           1        579.529000          1   579.529000
38           1        603.650000          1   603.650000
39           1        666.339000          1   666.339000
40           1        550.327000          1   550.327000
41           1        591.559000          1   591.559000
42           1        591.382000          1   591.382000
43           1        519.233000          1   519.233000
44           1        538.591000          1   538.591000
45           1        638.973000          1   638.973000
46           1        587.038000          1   587.038000
47           1        573.654000          1   573.654000
48           1        630.891000          1   630.891000
49           1        637.094000          1   637.094000
50           1        717.768000          1   717.768000
51           1        718.191000          1   718.191000
52           1        666.407000          1   666.407000
53           1        680.065000          1   680.065000
54           1        664.995000          1   664.995000
55           1        568.874000          1   568.874000
56           1        622.518000          1   622.518000
57           1        681.942000          1   681.942000
58           1        642.889000          1   642.889000
59           1        594.256000          1   594.256000
60           1        662.824000          1   662.824000
61           1        643.589000          1   643.589000
62           1        708.739000          1   708.739000
63           1        686.351000          1   686.351000
64           1        694.189000          1   694.189000
65           1        736.061000          1   736.061000
66           1        689.778000          1   689.778000
67           1        651.632000          1   651.632000
68           1        630.938000          1   630.938000
69           1        682.784000          1   682.784000
70           1        650.005000          1   650.005000
71           1        655.940000          1   655.940000
72           1        744.679000          1   744.679000
73           1        683.629000          1   683.629000
74           1        736.540000          1   736.540000
75           1        737.157000          1   737.157000
76           1        700.399000          1   700.399000
77           1        745.463000          1   745.463000
78           1        681.199000          1   681.199000
79           1        680.973000          1   680.973000
80           1        642.319000          1   642.319000
81           1        692.239000          1   692.239000
82           1        719.967000          1   719.967000
83           1        673.230000          1   673.230000
84           1        756.261000          1   756.261000
85           1        703.517000          1   703.517000
86           1        743.147000          1   743.147000
87           1        709.817000          1   709.817000
88           1        779.781000          1   779.781000
89           1        792.793000          1   792.793000
90           1        779.590000          1   779.590000
91           1        714.458000          1   714.458000
92           1        631.065000          1   631.065000
93           1        796.818000          1   796.818000
94           1        744.934000          1   744.934000
95           1        644.320000          1   644.320000
96           1        856.526000          1   856.526000
97           1        813.332000          1   813.332000
98           1        860.099000          1   860.099000
99           1        820.510000          1   820.510000
100          1        728.614000          1   728.614000
101          1        818.311000          1   818.311000
102          1        883.431000          1   883.431000
103          1        667.977000          1   667.977000
104          1        684.333000          1   684.333000
105          1        866.689000          1   866.689000
106          1        745.683000          1   745.683000
107          1        709.176000          1   709.176000
108          1        913.250000          1   913.250000
109          1        891.660000          1   891.660000
110          1        915.759000          1   915.759000
111          1        923.341000          1   923.341000
112          1        863.405000          1   863.405000
113          1        883.316000          1   883.316000
114          1        940.194000          1   940.194000
115          1        783.616000          1   783.616000
116          1        857.963000          1   857.963000
117          1        983.934000          1   983.934000
118          1        896.433000          1   896.433000
119          1        766.570000          1   766.570000
120          1        937.017000          1   937.017000
121          1        886.522000          1   886.522000
122          1       1000.110000          1  1000.110000
123          1        931.424000          1   931.424000
124          1        850.601000          1   850.601000
125          1        925.099000          1   925.099000
126          1        905.242000          1   905.242000
127          1        775.562000          1   775.562000
128          1        852.083000          1   852.083000
129          1        913.512000          1   913.512000
130          1        925.351000          1   925.351000
131          1        972.501000          1   972.501000
132          1        875.848000          1   875.848000
133          1        939.471000          1   939.471000
134          1       1025.930000          1  1025.930000
135          1        969.174000          1   969.174000
136          1        919.219000          1   919.219000
137          1       1025.520000          1  1025.520000
138          1        897.765000          1   897.765000
139          1        919.965000          1   919.965000
140          1        983.755000          1   983.755000
141          1        955.468000          1   955.468000
142          1       1017.850000          1  1017.850000
143          1        969.072000          1   969.072000
144          1       1031.050000          1  1031.050000
145          1       1084.380000          1  1084.380000
146          1       1127.130000          1  1127.130000
147          1       1003.800000          1  1003.800000
148          1       1140.560000          1  1140.560000
149          1        975.003000          1   975.003000
150          1       1084.030000          1  1084.030000
151          1       1032.700000          1  1032.700000
152          1        917.901000          1   917.901000
153          1       1195.890000          1  1195.890000
154          1       1113.430000          1  1113.430000
155          1        983.923000          1   983.923000
156          1       1183.240000          1  1183.240000
157          1       1190.710000          1  1190.710000
158          1       1228.540000          1  1228.540000
159          1       1388.620000          1  1388.620000
160          1       1167.270000          1  1167.270000
161          1       1130.490000          1  1130.490000
162          1       1280.810000          1  1280.810000
163          1       1039.620000          1  1039.620000
164          1       1059.450000          1  1059.450000
165          1       1271.790000          1  1271.790000
166          1       1157.880000          1  1157.880000
167          1        970.603000          1   970.603000
168          1       1302.330000          1  1302.330000
169          1       1198.840000          1  1198.840000
170          1       1129.150000          1  1129.150000
171          1       1344.040000          1  1344.040000
172          1       1174.650000          1  1174.650000
173          1       1200.870000          1  1200.870000
174          1       1276.140000          1  1276.140000
175          1        966.832000          1   966.832000
176          1       1163.050000          1  1163.050000
177          1       1319.850000          1  1319.850000
178          1       1117.500000          1  1117.500000
179          1        975.780000          1   975.780000
180          1       1293.250000          1  1293.250000
181          1       1229.270000          1  1229.270000
182          1       1256.380000          1  1256.380000
183          1       1290.970000          1  1290.970000
184          1       1184.200000          1  1184.200000
185          1       1235.310000          1  1235.310000
186          1       1210.460000          1  1210.460000
187          1       1026.420000          1  1026.420000
188          1       1175.950000          1  1175.950000
189          1       1360.120000          1  1360.120000
190          1       1216.240000          1  1216.240000
191          1       1111.180000          1  1111.180000
192          1       1279.310750          1  1279.310750
193          1       1331.310108          1  1331.310108
194          1       1431.796445          1  1431.796445
195          1       1416.135144          1  1416.135144
196          1       1264.227963          1  1264.227963
197          1       1329.251516          1  1329.251516
198          1       1242.510175          1  1242.510175
199          1       1189.446839          1  1189.446839
200          1       1324.179675          1  1324.179675
201          1       1420.271787          1  1420.271787
202          1       1378.022280          1  1378.022280
203          1       1257.538376          1  1257.538376
204          1       1578.826008          1  1578.826008
205          1       1423.930732          1  1423.930732
206          1       1515.827837          1  1515.827837
207          1       1604.634231          1  1604.634231
208          1       1486.693644          1  1486.693644
209          1       1569.316994          1  1569.316994
210          1       1418.075317          1  1418.075317
211          1       1284.135777          1  1284.135777
212          1       1426.166921          1  1426.166921
213          1       1626.197718          1  1626.197718
214          1       1573.925977          1  1573.925977
215          1       1304.485318          1  1304.485318
216          1       1693.664243          1  1693.664243
217          1       1584.197336          1  1584.197336
218          1       1758.798356          1  1758.798356
219          1       1591.510298          1  1591.510298
220          1       1592.936393          1  1592.936393
221          1       1518.670886          1  1518.670886
222          1       1555.355652          1  1555.355652
223          1       1316.781107          1  1316.781107
224          1       1384.224101          1  1384.224101
225          1       1640.951980          1  1640.951980
226          1       1630.405343          1  1630.405343
227          1       1154.201685          1  1154.201685
228          1       1802.369136          1  1802.369136
229          1       1573.444253          1  1573.444253
230          1       1572.894386          1  1572.894386
231          1       1726.117829          1  1726.117829
232          1       1584.012727          1  1584.012727
233          1       1541.813156          1  1541.813156
234          1       1605.865009          1  1605.865009
235          1       1347.606746          1  1347.606746
236          1       1549.537977          1  1549.537977
237          1       1742.708346          1  1742.708346
238          1       1675.316602          1  1675.316602
239          1       1383.203759          1  1383.203759
240          1       1934.114045          1  1934.114045
241          1       1847.766762          1  1847.766762
242          1       1997.396465          1  1997.396465
243          1       1864.006483          1  1864.006483
244          1       1649.612729          1  1649.612729
245          1       1877.115229          1  1877.115229
246          1       1702.329274          1  1702.329274
247          1       1457.675005          1  1457.675005
248          1       1764.000101          1  1764.000101
249          1       1788.872527          1  1788.872527
250          1       1855.210489          1  1855.210489
251          1       1519.451912          1  1519.451912
252          1       1745.458015          1  1745.458015
253          1       1948.216502          1  1948.216502
254          1       2147.177108          1  2147.177108
255          1       2042.870492          1  2042.870492
256          1       1799.633829          1  1799.633829
257          1       1986.671712          1  1986.671712
258          1       1868.842812          1  1868.842812
259          1       1630.148869          1  1630.148869
260          1       1790.836747          1  1790.836747
261          1       1965.865126          1  1965.865126
262          1       1814.941991          1  1814.941991
263          1       1607.625407          1  1607.625407
264          1       2056.250979          1  2056.250979
265          1       2026.678260          1  2026.678260
266          1       2316.554618          1  2316.554618
267          1       1997.301113          1  1997.301113
268          1       1895.008124          1  1895.008124
269          1       1915.037265          1  1915.037265
270          1       1848.730688          1  1848.730688
271          1       1703.809052          1  1703.809052
272          1       1958.588599          1  1958.588599
273          1       1935.128576          1  1935.128576
274          1       2050.561711          1  2050.561711
275          1       1788.185888          1  1788.185888
276          1       2093.944246          1  2093.944246
277          1       2010.011145          1  2010.011145
278          1       2190.960233          1  2190.960233
279          1       1821.988816          1  1821.988816
280          1       2160.473723          1  2160.473723
281          1       2161.584872          1  2161.584872
282          1       1923.413218          1  1923.413218
283          1       1872.786185          1  1872.786185
284          1       2066.577865          1  2066.577865
285          1       1984.327788          1  1984.327788
286          1       2079.146153          1  2079.146153
287          1       1677.632998          1  1677.632998
288          1       2218.696319          1  2218.696319
289          1       2165.900248          1  2165.900248
290          1       2278.931128          1  2278.931128
291          1       2125.671754          1  2125.671754
292          1       2269.975256          1  2269.975256
293          1       2216.538171          1  2216.538171
294          1       2341.314302          1  2341.314302
295          1       1845.753527          1  1845.753527
296          1       2102.092843          1  2102.092843
297          1       2274.554017          1  2274.554017
298          1       2276.984917          1  2276.984917
299          1       1837.492490          1  1837.492490
300          1       2564.425189          1  2564.425189
301          1       2316.882825          1  2316.882825
302          1       2319.543500          1  2319.543500
303          1       2836.134666          1  2836.134666
304          1       2298.994540          1  2298.994540
305          1       2628.123791          1  2628.123791
306          1       2601.888277          1  2601.888277
307          1       2011.319582          1  2011.319582
308          1       2467.638172          1  2467.638172
309          1       2800.181911          1  2800.181911
310          1       2518.641792          1  2518.641792
311          1       2283.650749          1  2283.650749
312          1       2503.325959          1  2503.325959
313          1       2780.621450          1  2780.621450
314          1       2955.455282          1  2955.455282
315          1       2840.789749          1  2840.789749
316          1       2497.314067          1  2497.314067
317          1       2742.423105          1  2742.423105
318          1       2758.012501          1  2758.012501
319          1       2061.011312          1  2061.011312
320          1       2679.369371          1  2679.369371
321          1       2673.709732          1  2673.709732
322          1       2597.637655          1  2597.637655
323          1       2212.575468          1  2212.575468
324          1       2432.166739          1  2432.166739
325          1       2604.593563          1  2604.593563
326          1       3233.568729          1  3233.568729
327          1       2699.130452          1  2699.130452
328          1       2662.002223          1  2662.002223
329          1       2884.175911          1  2884.175911
330          1       2768.304153          1  2768.304153
331          1       2381.557779          1  2381.557779
332          1       3098.683354          1  3098.683354
333          1       3221.334469          1  3221.334469
334          1       3245.327728          1  3245.327728
335          1       2759.946376          1  2759.946376
336          1       2871.436763          1  2871.436763
337          1       3268.689740          1  3268.689740
338          1       3396.638796          1  3396.638796
339          1       2568.870919          1  2568.870919
340          1       3368.824764          1  3368.824764
341          1       2917.814125          1  2917.814125
342          1       2894.924972          1  2894.924972
343          1       2863.185568          1  2863.185568
344          1       3076.542609          1  3076.542609
345          1       3063.501851          1  3063.501851
346          1       3271.496810          1  3271.496810
347          1       2329.749654          1  2329.749654
348          1       3620.523034          1  3620.523034
349          1       3472.109371          1  3472.109371
350          1       3789.560479          1  3789.560479
351          1       3452.247630          1  3452.247630
352          1       3878.813856          1  3878.813856
353          1       3576.276369          1  3576.276369
354          1       3631.762755          1  3631.762755
355          1       3144.688436          1  3144.688436
356          1       3223.053233          1  3223.053233
357          1       3871.772327          1  3871.772327
358          1       3402.581118          1  3402.581118
359          1       2769.154674          1  2769.154674
360          1       4177.105546          1  4177.105546
361          1       3596.429234          1  3596.429234
362          1       3728.505419          1  3728.505419
363          1       4247.070342          1  4247.070342
364          1       3930.391397          1  3930.391397
365          1       3901.663325          1  3901.663325
366          1       3942.881064          1  3942.881064
367          1       3408.622870          1  3408.622870
368          1       3726.807920          1  3726.807920
369          1       4017.764998          1  4017.764998
370          1       3410.331302          1  3410.331302
371          1       2856.014070          1  2856.014070
372          1       3752.336074          1  3752.336074
373          1       4088.540440          1  4088.540440
374          1       4060.577785          1  4060.577785
375          1       3596.428902          1  3596.428902
376          1       3802.155006          1  3802.155006
377          1       3908.611274          1  3908.611274
378          1       3947.581666          1  3947.581666
379          1       2976.529710          1  2976.529710
380          1       3893.696299          1  3893.696299
381          1       4323.885698          1  4323.885698
382          1       3732.355136          1  3732.355136
383          1       3110.941953          1  3110.941953
384          1       4282.664822          1  4282.664822
385          1       4311.081845          1  4311.081845
386          1       4569.473832          1  4569.473832
387          1       4131.874848          1  4131.874848
388          1       4082.986913          1  4082.986913
389          1       4340.404966          1  4340.404966
390          1       4106.358856          1  4106.358856
391          1       3622.355496          1  3622.355496
392          1       4331.565932          1  4331.565932
393          1       4269.316133          1  4269.316133
394          1       4357.228471          1  4357.228471
395          1       3196.612532          1  3196.612532
396          1       4552.491027          1  4552.491027
397          1       4404.788177          1  4404.788177
398          1       4691.859487          1  4691.859487
399          1       4466.672442          1  4466.672442
400          1       4657.794227          1  4657.794227
401          1       5009.121707          1  5009.121707
402          1       4374.383025          1  4374.383025
403          1       4181.309120          1  4181.309120
404          1       4870.324492          1  4870.324492
405          1       4777.218407          1  4777.218407
406          1       4830.374761          1  4830.374761
407          1       4021.645331          1  4021.645331
408          1       5484.585694          1  5484.585694
409          1       4880.564238          1  4880.564238
410          1       5636.832836          1  5636.832836
411          1       4822.765701          1  4822.765701
412          1       5406.392156          1  5406.392156
413          1       5346.499075          1  5346.499075
414          1       5268.496570          1  5268.496570
415          1       4894.786647          1  4894.786647
416          1       5216.915842          1  5216.915842
417          1       5934.567110          1  5934.567110
418          1       5703.432960          1  5703.432960
419          1       4379.032648          1  4379.032648
420          1       6377.001675          1  6377.001675
421          1       5525.170945          1  5525.170945
422          1       5904.197643          1  5904.197643
423          1       5419.860549          1  5419.860549
424          1       6213.531037          1  6213.531037
425          1       5837.081682          1  5837.081682
426          1       6203.451434          1  6203.451434
427          1       5185.086987          1  5185.086987
428          1       5658.695351          1  5658.695351
429          1       6112.263520          1  6112.263520
430          1       6167.840189          1  6167.840189
431          1       4206.727284          1  4206.727284
432          1       6010.403025          1  6010.403025
433          1       6277.912088          1  6277.912088
434          1       6065.605558          1  6065.605558
435          1       6679.369777          1  6679.369777
436          1       6196.796828          1  6196.796828
437          1       6562.105823          1  6562.105823
438          1       6993.545427          1  6993.545427
439          1       5303.806263          1  5303.806263
440          1       5853.143021          1  5853.143021
441          1       6151.706530          1  6151.706530
442          1       5444.528385          1  5444.528385
443          1       4379.356780          1  4379.356780
444          1       6157.785988          1  6157.785988
445          1       5597.823973          1  5597.823973
446          1       6013.368154          1  6013.368154
447          1       6190.533533          1  6190.533533
448          1       5835.120482          1  5835.120482
449          1       5767.473267          1  5767.473267
450          1       6690.408664          1  6690.408664
451          1       5484.044354          1  5484.044354
452          1       6062.546108          1  6062.546108
453          1       6495.015837          1  6495.015837
454          1       6323.374031          1  6323.374031
455          1       5153.749752          1  5153.749752
456          1       6581.277348          1  6581.277348
457          1       6011.557235          1  6011.557235
458          1       7322.960557          1  7322.960557
459          1       6586.340029          1  6586.340029
460          1       6422.818168          1  6422.818168
461          1       6473.321802          1  6473.321802
462          1       6692.573142          1  6692.573142
463          1       5791.215160          1  5791.215160
464          1       6000.860697          1  6000.860697
465          1       6135.724399          1  6135.724399
466          1       6529.411163          1  6529.411163
467          1       5361.333128          1  5361.333128
468          1       6327.800512          1  6327.800512
469          1       6396.310306          1  6396.310306
470          1       6963.410181          1  6963.410181
471          1       5821.055214          1  5821.055214
472          1       7490.113912          1  7490.113912
473          1       5601.896958          1  5601.896958

@jessegrabowski
Copy link
Owner

jessegrabowski commented Aug 7, 2024

So the problem is either in the covariance matrix or the loss function.

The covariance matrix should be easier to test, so I would start there. Can you check that for given rho, sigma, and n, build_chao_lin_covariance matches the R equivalent?

@steveshaoucsb
Copy link
Contributor Author

steveshaoucsb commented Aug 7, 2024

I only generated the chow-lin covariance from the python side as I am unable to obtain the one through the source code in R(looks unaccessible on my side). I've add that file as tests/data/chowlin_covariance.csv.
Also for your info(results from python package):
rho: 0.99999
sigma_e_sq: 4148.985003954791
n: 474

@jessegrabowski jessegrabowski self-requested a review August 8, 2024 01:50
@jessegrabowski
Copy link
Owner

I believe I fixed the bug. I misunderstood this line in the paper:

In a different approach, Bournay and Laroque (1979, p. 23) suggest the maximization of the likelihood of the GLS-regression:

$$ L(\rho, \sigma_{\epsilon}^2, \beta) = \frac{\exp \left[ -\frac{1}{2} u_1' (C \Sigma C')^{-1} u_1 \right]}{(2\pi)^{n_1/2} \cdot [\det (C \Sigma C')]^{1/2}}, $$

where $u_1$ is given by Eq. (2) and (6). $\hat{\beta}$ turns out to be the GLS estimator from Eq. (5). The maximum likelihood estimator of the autoregressive parameter, $\hat{\rho}$, is a consistent estimator of the true value, thus it has been chosen as the default estimator. However, in some cases, $\hat{\rho}$ turns out to be negative even if the true $\rho$ is positive. Thus, by default, tempdisagg constrains the optimization space for $\rho$ to positive values.

For some reason I took this to mean that we don't minimize beta as part of this objective function, and instead compute it using the values of rho and sigma obtained from the minimization. Looking back I have no idea how I reached this conclusion. But now it's fixed, and the outputs match R.

@jessegrabowski
Copy link
Owner

Thanks again for finding this, and for all your work debugging!

@jessegrabowski jessegrabowski merged commit 5d7d48f into jessegrabowski:main Aug 8, 2024
2 checks passed
@steveshaoucsb
Copy link
Contributor Author

No problem! Hopefully everything will work well after release!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants