1 /*
2 * Copyright (C) 1996-2011 Daniel Waggoner and Tao Zha
3 *
4 * This free software: you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation, either version 3 of the License, or
7 * (at your option) any later version.
8 *
9 * It is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * If you did not received a copy of the GNU General Public License
15 * with this software, see <http://www.gnu.org/licenses/>.
16 */
17
18 #include "dw_switch_opt.h"
19 #include "dw_std.h"
20
21 #include <stdlib.h>
22 #include <string.h>
23 #include <math.h>
24
25 //#define DW_NPSOL_PRESENT
26 //#define DW_CSMINWEL_PRESENT
27
28 //====== Static Global Variables ======
29 static TStateModel *Model=(TStateModel*)NULL;
30 static PRECISION *buffer=(PRECISION*)NULL;
31 static PRECISION *ModifiedFreeParameters=(PRECISION*)NULL;
32 static PRECISION *FreeParameters_Q=(PRECISION*)NULL;
33 static int NumberFreeParameters_Q=0;
34 static PRECISION *FreeParameters_Theta=(PRECISION*)NULL;
35 static int NumberFreeParameters_Theta=0;
36
37 /*
38 Sets up global variables for optimization. The buffer used to store the free
39 parameters is allocated and deallocated outside these routines. The buffers
40 FreeQ and FreeTheta do not have to be contiguous, but FreeQ must be at least
41 NumberFreeParametersQ(model)*sizeof(PRECISION) bytes and FreeTheta must be at
42 least NumberFreeParametersTheta(model)*sizeof(PRECISION) bytes. The region
43 from Modified + 0 through Modified + n - 1 must be completely contained in
44 FreeQ and FreeTheta. Here n is the parameter passed to MLEObjectiveFunction()
45 or PosteriorObjectiveFunction().
46 */
SetupObjectiveFunction(TStateModel * model,PRECISION * Modified,PRECISION * FreeQ,PRECISION * FreeTheta)47 void SetupObjectiveFunction(TStateModel *model, PRECISION *Modified, PRECISION *FreeQ, PRECISION *FreeTheta)
48 {
49 if (buffer) dw_free(buffer);
50 buffer=(PRECISION*)NULL;
51 Model=model;
52 FreeParameters_Q=FreeQ;
53 NumberFreeParameters_Q=NumberFreeParametersQ(model);
54 FreeParameters_Theta=FreeTheta;
55 NumberFreeParameters_Theta=NumberFreeParametersTheta(model);
56 ModifiedFreeParameters=Modified;
57 }
58
59 /*
60 Sets up global variables for optimization. The buffer used to store the free
61 parameters is allocated and deallocated in this routine. One of FreeTheta_Idx
62 or FreeQ_Idx must be zero and the other must be NumberFreeParametersQ(model)
63 or NumberFreeParameters(Theta). Modified_Idx and Modified_Idx + n -1 must be
64 between 0 and NumberFreeParametersQ(model) + NumberFreeParameters(Theta) - 1,
65 inclusive. Here n is the parameter passed to PosteriorObjectiveFunction() or
66 MLEObjectiveFunction().
67
68 */
SetupObjectiveFunction_AllocateMemory(TStateModel * model,int FreeTheta_Idx,int FreeQ_Idx,int Modified_Idx)69 void SetupObjectiveFunction_AllocateMemory(TStateModel *model, int FreeTheta_Idx, int FreeQ_Idx, int Modified_Idx)
70 {
71 if (buffer) dw_free(buffer);
72 Model=model;
73 NumberFreeParameters_Q=NumberFreeParametersQ(model);
74 NumberFreeParameters_Theta=NumberFreeParametersTheta(model);
75 buffer=(PRECISION*)dw_malloc((NumberFreeParameters_Q + NumberFreeParameters_Theta)*sizeof(PRECISION));
76
77 FreeParameters_Q=buffer+FreeQ_Idx;
78 FreeParameters_Theta=buffer+FreeTheta_Idx;
79 ModifiedFreeParameters=buffer+Modified_Idx;
80 }
81
PosteriorObjectiveFunction(PRECISION * x,int n)82 PRECISION PosteriorObjectiveFunction(PRECISION *x, int n)
83 {
84 if (x != ModifiedFreeParameters) memmove(ModifiedFreeParameters,x,n*sizeof(PRECISION));
85 ConvertFreeParametersToQ(Model,FreeParameters_Q);
86 ConvertFreeParametersToTheta(Model,FreeParameters_Theta);
87 return -LogPosterior_StatesIntegratedOut(Model);
88 }
89
PosteriorObjectiveFunction_csminwel(double * x,int n,double ** args,int * dims)90 PRECISION PosteriorObjectiveFunction_csminwel(double *x, int n, double **args, int *dims)
91 {
92 return PosteriorObjectiveFunction(x,n);
93 }
94
PosteriorObjectiveFunction_npsol(int * mode,int * n,double * x,double * f,double * g,int * nstate)95 void PosteriorObjectiveFunction_npsol(int *mode, int *n, double *x, double *f, double *g, int *nstate)
96 {
97 *f=PosteriorObjectiveFunction(x,*n);
98 }
99
MLEObjectiveFunction(PRECISION * x,int n)100 PRECISION MLEObjectiveFunction(PRECISION *x, int n)
101 {
102 if (x != ModifiedFreeParameters) memmove(ModifiedFreeParameters,x,n*sizeof(PRECISION));
103 ConvertFreeParametersToQ(Model,FreeParameters_Q);
104 ConvertFreeParametersToTheta(Model,FreeParameters_Theta);
105 return -LogLikelihood_StatesIntegratedOut(Model);
106 }
107
MLEObjectiveFunction_csminwel(double * x,int n,double ** args,int * dims)108 PRECISION MLEObjectiveFunction_csminwel(double *x, int n, double **args, int *dims)
109 {
110 return MLEObjectiveFunction(x,n);
111 }
112
MLEObjectiveFunction_npsol(int * mode,int * n,double * x,double * f,double * g,int * nstate)113 void MLEObjectiveFunction_npsol(int *mode, int *n, double *x, double *f, double *g, int *nstate)
114 {
115 *f=MLEObjectiveFunction(x,*n);
116 }
117
118 #ifdef DW_NPSOL_PRESENT
119 //=== NPSOL declarations ===
120 #define npsol npsol_
121 #define npoptn npoptn_
122 #define npsol_options npsol_options_
123 void npsol(int *n, int *nclin, int *ncnln, int *ldA, int *ldJ, int *ldR, double *A, double *bl, double *bu,
124 void (*funcon)(int *mode, int *ncnln, int *n, int *ldJ, int *needc, double *x, double *c, double *cJac, int *nstate),
125 void (*funobj)(int *mode, int *n, double *x, double *f, double *g, int *nstate),
126 int *inform, int *iter, int *istate, double *c, double *cJac, double *clamda, double *f, double *g, double *R, double *x,
127 int *iw, int *leniw, double *w, int *lenw);
128 void npoptn(char *,int);
confun(int * mode,int * ncnln,int * n,int * ldJ,int * needc,double * x,double * c,double * cJac,int * nstate)129 static void confun(int *mode, int *ncnln, int *n, int *ldJ, int *needc, double *x, double *c, double *cJac, int *nstate)
130 {
131 *mode=-1;
132 }
PosteriorOptimizationQ_npsol(TStateModel * model,int max_iteration)133 int PosteriorOptimizationQ_npsol(TStateModel *model, int max_iteration)
134 {
135 int j, size_theta, pos_theta, size_Q, pos_Q, opt;
136 double objective, objective_last, likelihood, prior;
137 char option_string[256];
138
139 // npsol arguments
140 int n, nclin=0, ncnln=0, ldA=1, ldJ=1, ldR, nctotl, *istate, *iw, leniw, lenw, inform, iter;
141 double A[1], cJac[1], *bl, *bu, *clamda, *R, *x, *w, f, *g, c[1], tolerance;
142
143 //==== Allocate memory ===
144 size_theta=NumberFreeParametersTheta(model);
145 size_Q=NumberFreeParametersQ(model);
146 pos_theta=0;
147 pos_Q=size_theta;
148
149 n=size_Q;
150 nctotl=n+nclin+ncnln;
151 leniw=3*n + nclin + 2*ncnln;
152 if (ncnln == 0)
153 if (nclin == 0)
154 lenw=20*n;
155 else
156 lenw=2*n*n + 11*nclin;
157 else
158 lenw=2*n*n + n*nclin + 2*n*ncnln + 20*n + 11*nclin + 21*ncnln;
159 istate=(int*)dw_malloc(nctotl*sizeof(int));
160 iw=(int*)dw_malloc(leniw*sizeof(int));
161 bl=(double*)dw_malloc(nctotl*sizeof(double));
162 bu=(double*)dw_malloc(nctotl*sizeof(double));
163 clamda=(double*)dw_malloc(nctotl*sizeof(double));
164 R=(double*)dw_malloc(n*n*sizeof(double));
165 x=(double*)dw_malloc((size_Q + size_theta)*sizeof(double));
166 g=(double*)dw_malloc(n*sizeof(double));
167 w=(double*)dw_malloc(lenw*sizeof(double));
168 //====================================================================
169
170 //=== Set starting value ===
171 ConvertQToFreeParameters(model,x+pos_Q);
172 ConvertThetaToFreeParameters(model,x+pos_theta);
173
174 SetupObjectiveFunction(model,x+pos_Q,x+pos_Q,x+pos_theta);
175 ldR=n;
176 for (j=0; j < n; j++)
177 {
178 bl[j]=0.0;
179 bu[j]=1.0;
180 }
181
182 //=== Set NPSOL options
183 npoptn("Derivative level = 0",20);
184 sprintf(option_string,"Optimality tolerance = %le",tolerance);
185 npoptn(option_string,strlen(option_string));
186 sprintf(option_string,"Major iterations limit = %d",max_iteration);
187 npoptn(option_string,strlen(option_string));
188
189 npsol(&n,&nclin,&ncnln,&ldA,&ldJ,&ldR,A,bl,bu,&confun,&PosteriorObjectiveFunction_npsol,
190 &inform,&iter,istate,c,cJac,clamda,&f,g,R,x+pos_Q,iw,&leniw,w,&lenw);
191
192 ConvertFreeParametersToQ(model,x+pos_Q);
193 ConvertFreeParametersToTheta(model,x+pos_theta);
194
195 //=== Free memory ===
196 dw_free(istate);
197 dw_free(iw);
198 dw_free(bl);
199 dw_free(bu);
200 dw_free(clamda);
201 dw_free(R);
202 dw_free(x);
203 dw_free(g);
204 dw_free(w);
205 }
206 #endif
207
208 /* #ifdef DW_CSMINWEL_PRESENT */
209 /* /\* */
210 /* Maximizes the posterior kernel with respect to a subblock of theta or q given */
211 /* by offset and length. It must be the case that offset plus length be less than */
212 /* or equal to the number of free parameters in theta or q. If flag is 1, then */
213 /* subblocks of theta are used; if flag is 2, then subblocks of q are used, if */
214 /* flag is 3, then maximization over both theta and q is performed. */
215
216 /* Returns 1 upon completion and 0 if the arguments are invalid. */
217 /* *\/ */
218 /* int PosteriorOptimization_csminwel(TStateModel *model, int offset, int length, int flag, int max_iteration, PRECISION tolerance) */
219 /* { */
220 /* int size_theta, pos_theta, size_Q, pos_Q; */
221
222 /* // csminwel arguments */
223 /* int itct, fcount, retcodeh; */
224 /* double *x, fh; */
225 /* TMatrix H; */
226 /* TVector g; */
227
228 /* //==== Sizes */
229 /* size_theta=NumberFreeParametersTheta(model); */
230 /* size_Q=NumberFreeParametersQ(model); */
231 /* pos_theta=0; */
232 /* pos_Q=size_theta; */
233 /* switch (flag) */
234 /* { */
235 /* case 1: */
236 /* offset+=pos_theta; */
237 /* if (offset+length > size_theta) return 0; */
238 /* break; */
239 /* case 2: */
240 /* offset+=pos_Q; */
241 /* if (offset+length > size_Q) return 0; */
242 /* break; */
243 /* case 3: */
244 /* offset=0; */
245 /* length=size_Q+size_theta; */
246 /* } */
247
248 /* //=== Allocate memory */
249 /* x=(PRECISION*)dw_malloc((size_Q + size_theta)*sizeof(double)); */
250 /* g=CreateVector(length); */
251 /* H=IdentityMatrix((TMatrix)NULL,length); */
252 /* ProductMS(H,H,1e-5); */
253
254 /* //=== Setup objective function */
255 /* SetupObjectiveFunction(model,x+offset,x+pos_Q,x+pos_theta); */
256
257 /* //=== Set starting value */
258 /* ConvertQToFreeParameters(model,x+pos_Q); */
259 /* ConvertThetaToFreeParameters(model,x+pos_theta); */
260
261 /* //=== Optimize */
262 /* csminwel(PosteriorObjectiveFunction_csminwel,x+offset,length,pElementM(H),pElementV(g),NULL, */
263 /* &fh,tolerance,&itct,max_iteration,&fcount,&retcodeh,NULL,NULL); */
264
265 /* //=== Get terminal values */
266 /* ConvertFreeParametersToQ(model,x+pos_Q); */
267 /* ConvertFreeParametersToTheta(model,x+pos_theta); */
268
269 /* //=== Free memory === */
270 /* FreeVector(g); */
271 /* FreeMatrix(H); */
272 /* dw_free(x); */
273
274 /* return 1; */
275 /* } */
276 /* #endif */
277
278 #ifdef DW_CSMINWEL_PRESENT
279
280 #include "csminwel.h"
281
282 /*
283 Maximizes the posterior kernel with respect to a subblock of theta or q given
284 by offset and length. It must be the case that offset plus length be less than
285 or equal to the number of free parameters in theta or q. If flag is 1, then
286 subblocks of theta are used; if flag is 2, then subblocks of q are used, if
287 flag is 3, then maximization over both theta and q is performed.
288
289 Returns 1 upon completion and 0 if the arguments are invalid.
290 */
PosteriorOptimization_csminwel(TStateModel * model,int offset,int length,int flag,int max_iteration,PRECISION tolerance,PRECISION stepsize,PRECISION hessian_scale)291 int PosteriorOptimization_csminwel(TStateModel *model, int offset, int length, int flag, int max_iteration,
292 PRECISION tolerance, PRECISION stepsize, PRECISION hessian_scale)
293 {
294 int size_theta, pos_theta, size_Q, pos_Q;
295 PRECISION *parameters;
296
297 // csminwel arguments
298 int itct, fcount, retcodeh;
299 double *x, fh, original_stepsize=GRADSTPS_CSMINWEL;
300 TMatrix H;
301 TVector g;
302
303 //==== Sizes
304 size_theta=NumberFreeParametersTheta(model);
305 size_Q=NumberFreeParametersQ(model);
306 pos_theta=0;
307 pos_Q=size_theta;
308 if (flag == 3)
309 if (length > 0)
310 {
311 if (offset+length > size_theta + size_Q) return 0;
312 }
313 else
314 {
315 offset=0;
316 length=size_Q + size_theta;
317 }
318 else if (flag == 2)
319 if (length > 0)
320 {
321 offset+=pos_Q;
322 if (offset+length > size_Q) return 0;
323 }
324 else
325 {
326 offset=pos_Q;
327 length=size_Q;
328 }
329 else if (flag == 1)
330 if (length > 0)
331 {
332 offset+=pos_theta;
333 if (offset+length > size_theta) return 0;
334 }
335 else
336 {
337 offset=pos_theta;
338 length=size_theta;
339 }
340 else
341 {
342 printf("PosteriorOptimization_csminwel(): Unknown flag\n");
343 dw_exit(0);
344 }
345
346 //=== Allocate memory and set Hessian scale and step size
347 parameters=(PRECISION*)dw_malloc((size_Q + size_theta)*sizeof(double));
348 x=(PRECISION*)dw_malloc(length*sizeof(PRECISION));
349 g=CreateVector(length);
350 H=IdentityMatrix((TMatrix)NULL,length);
351 ProductMS(H,H,hessian_scale);
352 GRADSTPS_CSMINWEL=stepsize;
353
354 //=== Setup objective function
355 SetupObjectiveFunction(model,parameters+offset,parameters+pos_Q,parameters+pos_theta);
356
357 //=== Set starting value
358 ConvertQToFreeParameters(model,parameters+pos_Q);
359 ConvertThetaToFreeParameters(model,parameters+pos_theta);
360 memcpy(x,parameters+offset,length*sizeof(PRECISION));
361
362 //=== Optimize
363 csminwel(PosteriorObjectiveFunction_csminwel,x,length,pElementM(H),pElementV(g),NULL,
364 &fh,tolerance,&itct,max_iteration,&fcount,&retcodeh,NULL,NULL);
365
366 //=== Get terminal values
367 ConvertFreeParametersToQ(model,parameters+pos_Q);
368 ConvertFreeParametersToTheta(model,parameters+pos_theta);
369
370 //=== Free memory
371 FreeMatrix(H);
372 FreeVector(g);
373 dw_free(x);
374 dw_free(parameters);
375
376 //=== Reset global step size
377 GRADSTPS_CSMINWEL=original_stepsize;
378
379 return 1;
380 }
381
Estimate_csminwel(TStateModel * model,int iter_block,int iter_start,PRECISION iter_inc,PRECISION crit_start,PRECISION crit_end,PRECISION crit_inc,PRECISION stepsize_theta,PRECISION hessian_scale_theta,PRECISION stepsize_q,PRECISION hessian_scale_q,PRECISION stepsize_full,PRECISION hessian_scale_full,FILE * f_out)382 int Estimate_csminwel(TStateModel *model, int iter_block, int iter_start, PRECISION iter_inc,
383 PRECISION crit_start, PRECISION crit_end, PRECISION crit_inc,
384 PRECISION stepsize_theta, PRECISION hessian_scale_theta,
385 PRECISION stepsize_q, PRECISION hessian_scale_q,
386 PRECISION stepsize_full, PRECISION hessian_scale_full,FILE *f_out)
387 {
388 int nq, ntheta, i, total_iteration=1, iteration, iter=iter_start;
389 PRECISION crit=crit_start, posterior, posterior_last, likelihood, likelihood_last, original_stepsize=GRADSTPS_CSMINWEL;
390 TVector theta=(TVector)NULL, q=(TVector)NULL;
391 time_t *currenttime;
392
393 nq=NumberFreeParametersQ(model);
394 ntheta=NumberFreeParametersTheta(model);
395 if (nq > 0) q=CreateVector(nq);
396 if (ntheta > 0) theta=CreateVector(ntheta);
397
398 posterior=LogPosterior_StatesIntegratedOut(model);
399
400 if (f_out)
401 {
402 fprintf(f_out,"\n\n//=== Initial Value ===//\n");
403 fprintf(f_out,"likelihood value: %22.14le\n",LogLikelihood_StatesIntegratedOut(model));
404 fprintf(f_out,"posterior value: %22.14le\n",posterior);
405 if (ntheta > 0)
406 {
407 fprintf(f_out,"theta:\n");
408 ConvertThetaToFreeParameters(model,pElementV(theta));
409 dw_PrintVector(f_out,theta,"%.13le ");
410 }
411 if (nq > 0)
412 {
413 fprintf(f_out,"q:\n");
414 ConvertQToFreeParameters(model,pElementV(q));
415 dw_PrintVector(f_out,q,"%.13le ");
416 }
417 fflush(f_out);
418 }
419
420 currenttime=(time_t*)dw_malloc(4*sizeof(time_t));
421
422 for ( ; crit >= crit_end; crit*=crit_inc, iter=(int)((PRECISION)iter*iter_inc))
423 {
424 for (iteration=1; iteration <= iter_block; iteration++)
425 {
426 posterior_last=posterior;
427 time(currenttime);
428
429 // optimize over theta for model i
430 PosteriorOptimization_csminwel(model,0,0,1,iter,crit,stepsize_theta,hessian_scale_theta);
431 time(currenttime+1);
432
433 // optimize over q
434 PosteriorOptimization_csminwel(model,0,0,2,iter,crit,stepsize_q,hessian_scale_q);
435 time(currenttime+2);
436
437 // full optimization
438 PosteriorOptimization_csminwel(model,0,0,3,iter,crit,stepsize_full,hessian_scale_full);
439 time(currenttime+3);
440
441 posterior=LogPosterior_StatesIntegratedOut(model);
442
443 if (f_out)
444 {
445 fprintf(f_out,"\n//=== Iteration %d ===//\n",total_iteration);
446 fprintf(f_out,"Criterion/Max Iteration: %le %d\n",crit,iter);
447 fprintf(f_out,"Elapsed time for theta: %0.4f\n",difftime(currenttime[1],currenttime[0]));
448 fprintf(f_out,"Elapsed time for q: %0.4f\n",difftime(currenttime[2],currenttime[1]));
449 fprintf(f_out,"Elapsed time for full: %0.4f\n",difftime(currenttime[3],currenttime[2]));
450 fprintf(f_out,"likelihood value: %22.14le\n",LogLikelihood_StatesIntegratedOut(model));
451 fprintf(f_out,"posterior value: %22.14le\n",posterior);
452 fprintf(f_out,"%s",ctime(¤ttime[3]));
453 if (ntheta > 0)
454 {
455 fprintf(f_out,"theta:\n");
456 ConvertThetaToFreeParameters(model,pElementV(theta));
457 dw_PrintVector(f_out,theta,"%.13le ");
458 }
459 if (nq > 0)
460 {
461 fprintf(f_out,"q:\n");
462 ConvertQToFreeParameters(model,pElementV(q));
463 dw_PrintVector(f_out,q,"%.13le ");
464 }
465 fflush(f_out);
466 }
467
468 total_iteration++;
469 if (fabs(posterior - posterior_last) <= crit) break;
470 }
471 }
472
473 FreeVector(q);
474 FreeVector(theta);
475 dw_free(currenttime);
476
477 return 1;
478 }
479 #endif
480
481
482
483