C++

Forum

 
Program received signal SIGSEGV, Segmentation fault

Pages: 123
lifeisgood (51)
Hi, I have encountered the above error while debugging. I am suspecting that I have used some pointers to matrices incorrectly, so I intend to check the code over again, but it would be great if I could get some expert opinion, since I am still pretty new to the program.

The code is an attempt to replicate a paper that uses value function iterations to solve for a decision rule.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
 
/*
* File: main.cpp
* Author: Lia Yin
*
* Created on June 5, 2017
*/

#include <iostream>
#include <vector>
#include <algorithm> // std::count
#include <functional>
#include <tuple>
#include <cmath>

/* Includes, system */
#include <iomanip>
#include <fstream>

// define the variables and parameters
#define nk 60
#define nb 80
#define ne 8
#define na 2
#define nr 2
#define np 2

int skip;

/*
 * used in the section for testing the convergence of decision rules
 */
bool IsNonzero(int i){
    return(i!=0);
} // if the number is not equal to 0, IsNonzero returns a 1

int main(){
    double test2; // criterion for the difference between V and VA
    double test2_a;

    int i; // this one and the ones below are counters
    int j;
    int l;
    int m;
    int n;
    int o;

    // these variables are for the VFs
    double tempa; // the value of TFP for a specific combination of shocks
    double tempr; // the value of world interest rate
    double tempqb; // the value of bond price
    double tempp; // the value of int'l goods price
    double tempk; // the value of capital
    double templ; // the value of labor
    double tempv; // the value of imported goods
    double tempy; // the value of output
    double tempi; // the value of investment

    double tempq; // the value of the price of equity
    double tempcolat; // the collateral constraint

    double templb; // labor with collateral constraint binding
    double tempvb; // foreign inputs under collateral constraint
    double tempyb; // output under collateral constraint

    double p1nval; // value in the parentheses of the utility function
    double p1n[nk][nb]={-999999999}; // for each pair of k' and b', what is the value function equal to

    double expv; // expected value

    int kkn[nk][nb][ne]; // counter for next period's capital
    int bbn[nk][nb][ne]; // counter for next period's bond
    int kzt[nk][nb][ne]; // counter for last period's capital
    int kz2[nk][nb][ne]; // counter for capital from two periods ago
    int bzt[nk][nb][ne]; // counter for last period's bond
    int bz2[nk][nb][ne]; // counter for bond from two periods ago

    int con5; // counter for breaking the iterations down into smaller pieces: in total there are 2000 iterations; break them down into pieces of 26 iterations each
    int con; // counter for total number of iterations; the number is 2000

    int jjkk[nk][nb][ne]; // tally up to see if this group of iterations has the same k' value as last group of iterations and the group of iterations from two times ago
    int jjbb[nk][nb][ne];

    int jjk=0;
    int jjb=0;

    bool binding; // whether the collateral constraint is binding

    double VA[nk][nb][ne]; // updated value after each value function iteration

    bool leave=false; // break out of the loop if necessary

/*
 * INITIAL VALUES FOR PARAMETERS
 */
    const double maxit=2000; // max VF iterations with first price guess
    const double full=2; // VF iterations maximizing over k,b grids out of a total of VF iterations defined in "ALL"
    const double all=26; // "all-full" VF iterations do not maximize, just iterate using VF & dec rules of previous iteration
    double switching=0.00000005; // sets value of "VFDIF" at which VF ITERATIONS SWITCH TO MAXIMIZATION AT ALL TIMES

/* Mexican calibration values */
    const double sigma=2; // CRRA coefficient
    const double rbar=0.08571; // world real interest rate
    const double omega=1.8461; // labor elasticity in the utility function
    const double alpha=0.5924043; // share of capital in the production function
    const double beta=0.3053667;
    const double eta=1-alpha-beta; // share of foreign goods in the production function
    const double phi=0.2579; // fraction of inputs of V and m financed by working capital from firm-level data
    const double A=6.982; // productivity coefficient
    const double kappa=1e-150; // ceiling on the leverage ratio; 1*10^(-149.5701605)
    const double a=2.75; // capital adjustment cost coefficient
    const double delta=0.08800443; // capital depreciation rate
    const double gamma=0.01660445; // semielasticity of the rate of time preference with respect to composite good c-N(L)

/* steady state values */
    const double abar=1; // steady state technology
    const double taxbar=0.169; // steady state tax
    const double pbar=1.028; // steady state price

    const double kbar = 799.592; // steady state capital
    const double lbar = 19.050; // steady state labor
    const double vbar = 45.241; // steady state foreign goods
    const double cbar=265.43703344; // steady state consumption

    double wel=cbar-pow(lbar,omega)/omega; // one component of the steady state utility function
    const double rate=exp(gamma*log(1+wel)); // steady state time preference
    const double wel1=pow(wel,1-sigma)/(1-sigma); // steady state period utility function
    wel=wel1*(1/(1-1/rate)); // steady state total expected utility function

    double V[nk][nb][ne]={wel}; // initial guess of the value function

/* specification of markov process of exogenous shocks */
    double B[nb]; // bond grid with nb elements
    double K[nk]; // capital grid with nk elements
    double E[ne]; // shock grid with ne elements
    double P[ne][ne]={0}; // matrix of transition probabilities, initial values of which are all 0.0
    double Ashock[2]; // vector of TFP shocks
    Ashock[0] = log(1-0.0134);Ashock[1] = log(1+0.0134);
    double Rshock[2]; // vector of int'l rate shocks
    Rshock[0] = log(1-0.0196);Rshock[1] = log(1+0.0196);
    double Pshock[2]; // vector of int'l goods price shocks
    Pshock[0] = log(1-0.0335);Pshock[1] = log(1+0.0335);
lifeisgood (51)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
 
/*
 * CALCULATE THE STATE TRANSITION MATRIX
 */

/* A. shocks in int'l goods prices are independent of shocks in TFP or shocks in int'l interest rate, so let's take care of this part first */
    double rhop = 0.737; // autocorrelation coefficient for int'l goods prices
    double pp = 0.5; // if we just consider P, what's the probability its shock is of a certain form
    double pip[2][2]; // transition probabilities between the two states of P
    pip[1][1]=(1-rhop)*pp + rhop;
    pip[2][1]=1-pip[1][1];
    pip[2][2]=(1-rhop)*pp + rhop;
    pip[1][2]=1-pip[2][2];

/* B. joint Markov process for TFP and R shocks */
    double rhoa = 0.537; // autocorrelation coefficient for TFP
    double rhor = 0.572; // autocorrelation coefficient for int'l rate
    double rhoar = (rhoa+rhor)/2; // the procedure requires that the AR(1) coefficients of the shocks that are correlated with each other be the same
    double corrar = -0.67; // correlation coefficient between TFP and int'l rate

    int nar = na*nr; // # of combinations between TFP and R
    double par[nar]; // the probabilities of the four different combinations of TFP and R
    par[1] = (1+corrar)/4; // the probability that both shocks are positive
    par[2]=0.5-par[1];
    par[3]=par[2];
    par[4]=par[1];

    double piar[4][4]; // transition probabilities for TFP and R combined
    int b; // a boolean
    for(i=0;i<4;i=i+1) {
        for(j=0;j<4;j=j+1) {
            if(i=j) {
                b=0;
            }
            else{
                b=1;
            }
            piar[i][j]=(1-rhoar)*par[i]+rhoar*b;
        }
    }

/* C. Combine all three things, TFP, R, and P, make the large transition probability matrix */
    double pi[8][8];
    pi[1][1]=pip[1][1]*piar[1][1];
    pi[2][1]=pip[1][1]*piar[2][1];
    pi[3][1]=pip[2][1]*piar[1][1];
    pi[4][1]=pip[2][1]*piar[2][1];
    pi[5][1]=pip[1][1]*piar[3][1];
    pi[6][1]=pip[1][1]*piar[4][1];
    pi[7][1]=pip[2][1]*piar[3][1];
    pi[8][1]=pip[2][1]*piar[4][1];
    pi[1][2]=pip[1][1]*piar[1][2];
    pi[2][2]=pip[1][1]*piar[2][2];
    pi[3][2]=pip[2][1]*piar[1][2];
    pi[4][2]=pip[2][1]*piar[2][2];
    pi[5][2]=pip[1][1]*piar[3][2];
    pi[6][2]=pip[1][1]*piar[4][2];
    pi[7][2]=pip[2][1]*piar[3][2];
    pi[8][2]=pip[2][1]*piar[4][2];
    pi[1][3]=pip[1][2]*piar[1][1];
    pi[2][3]=pip[1][2]*piar[2][1];
    pi[3][3]=pip[2][2]*piar[1][1];
    pi[4][3]=pip[2][2]*piar[2][1];
    pi[5][3]=pip[1][2]*piar[3][1];
    pi[6][3]=pip[1][2]*piar[4][1];
    pi[7][3]=pip[2][2]*piar[3][1];
    pi[8][3]=pip[2][2]*piar[4][1];
    pi[1][4]=pip[1][2]*piar[1][2];
    pi[2][4]=pip[1][2]*piar[2][2];
    pi[3][4]=pip[2][2]*piar[1][2];
    pi[4][4]=pip[2][2]*piar[2][2];
    pi[5][4]=pip[1][2]*piar[3][2];
    pi[6][4]=pip[1][2]*piar[4][2];
    pi[7][4]=pip[2][2]*piar[3][2];
    pi[8][4]=pip[2][2]*piar[4][2];

    pi[1][5]=pip[1][1]*piar[1][3];
    pi[2][5]=pip[1][1]*piar[2][2];
    pi[3][5]=pip[2][1]*piar[1][3];
    pi[4][5]=pip[2][1]*piar[2][3];
    pi[5][5]=pip[1][1]*piar[3][3];
    pi[6][5]=pip[1][1]*piar[4][3];
    pi[7][5]=pip[2][1]*piar[3][3];
    pi[8][5]=pip[2][1]*piar[4][3];
    pi[1][6]=pip[1][1]*piar[1][4];
    pi[2][6]=pip[1][1]*piar[2][4];
    pi[3][6]=pip[2][1]*piar[1][4];
    pi[4][6]=pip[2][1]*piar[2][4];
    pi[5][6]=pip[1][1]*piar[3][4];
    pi[6][6]=pip[1][1]*piar[4][4];
    pi[7][6]=pip[2][1]*piar[3][4];
    pi[8][6]=pip[2][1]*piar[4][4];
    pi[1][7]=pip[1][2]*piar[1][3];
    pi[2][7]=pip[1][2]*piar[2][3];
    pi[3][7]=pip[2][2]*piar[1][3];
    pi[4][7]=pip[2][2]*piar[2][3];
    pi[5][7]=pip[1][2]*piar[3][3];
    pi[6][7]=pip[1][2]*piar[4][3];
    pi[7][7]=pip[2][2]*piar[3][3];
    pi[8][7]=pip[2][2]*piar[4][3];
    pi[1][8]=pip[1][2]*piar[1][4];
    pi[2][8]=pip[1][2]*piar[2][4];
    pi[3][8]=pip[2][2]*piar[1][4];
    pi[4][8]=pip[2][2]*piar[2][4];
    pi[5][8]=pip[1][2]*piar[3][4];
    pi[6][8]=pip[1][2]*piar[4][4];
    pi[7][8]=pip[2][2]*piar[3][4];
    pi[8][8]=pip[2][2]*piar[4][4];

/* joint limiting distribution of R,P,TFP shocks */
    double temppi[8][8];
    for(i=0;i<8;i=i+1){
        for(j=0;j<8;j=j+1){
            temppi[i][j]=pi[i][j];
        }
    } // initialize temppi=pi

    double temppi_a[8][8];
    for(i=0;i<300;i=i+1){
        for(j=0;j<8;j=j+1){
            for(m=0;m<8;m=m+1){
                for(l=0;l<8;l=l+1){
                    temppi_a[j][m]=temppi[j][l]*pi[l][m]+temppi_a[j][m];
                    temppi[j][m]=temppi_a[j][m];
                }
            }
        } // finally finished matrix multiplication
        for(j=0;j<8;j=j+1){
            for(l=0;l<8;l=l+1){
                temppi[j][l]=temppi_a[j][l];
            }
        } // change the values of temppi for the next round of multiplication
    } // transition matrix powered 300 times
    double prob[ne];
    for(i=0;i<ne;i=i+1){
        prob[i]=temppi[i][1]; // becomes a row vector
    } // joint limiting probability of shocks
    m=0;
    double sa[8];
    double sr[8];
    double qb[8];
    double sp[8];
    for(i=0;i<nr;i=i+1){
        for(j=0;j<np;j=j+1){
            for(l=0;l<na;l=l+1){
                m=m+1;
                sa[m]=exp(Ashock[l])*abar; // state for TFP
                sr[m]=exp(Rshock[i])*(1+rbar); // state for the interest rate
                qb[m]=1/sr[m]; // state for bond prices
                sp[m]=exp(Pshock[j])*pbar; // state for price of int'l goods
            }
        }
    }

/*
 * CREATE GRIDS
 */

/* A. domestic capital */
    double kss=752.270;
    double kgridgap=3.2;
    double klow=kss-kgridgap*(nk/2);
    double k[nk];
    for(i=0;i<nk;i=i+1){
        k[i]=klow+i*kgridgap;
    }

/* B. foreign bonds */
    double bss=-127.184;
    double bgridgap=2*kgridgap;
    double blow=bss-bgridgap*(nb/2);
    double bond[nb];
    for(i=0;i<nb;i=i+1){
        bond[i]=blow+i*bgridgap;
    }
lifeisgood (51)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
 
/*
 * CALCULATE THE VALUE OF CHOOSING CAPITAL K'(M) AND ASSET HOLDINGS B'(N) WHEN THE INITIAL STATE IS
 * GIVEN BY (L,J,I)
 */

    skip=0;
    test2=switching;

/* A. calculate labor, int'l inputs allocations for given K(I), S(IS) assuming credit constraint does
 not bind; note that they are independent of K' and B' */

    for(i=0;i<ne;i=i+1){
        for(j=0;j<nb;j=j+1){
            for(l=0;l<nk;l=l+1){
                kzt[i][j][l]=-999999999;
                kz2[i][j][l]=-999999999;
            }
        }
    } // initialize the values of kzt and kz2 to make sure that we don't stop at the first iteration

    for(con5=0;con5<all;con5=con5+1){ // every 26 iterations
        con=con+1; // add one every time it loops
        while(con<maxit){ // whenever con<2000
            if(test2<switching){
                skip=1; // when V is similar to VA, restart maximization
            }
            for(i=0;i<ne;i=i+1){ // this period's shocks
                tempa=sa[i]; // the value of TFP for this specific combination of shocks
                tempr=sr[i]-1; // the value of world interest rate
                tempqb=qb[i]; // the value of bond price
                tempp=sp[i]; // the value of int'l goods price
                for(j=0;j<nb;j=j+1){ // this period's bond
                    for(l=0;l<nk;l=l+1){ // this period's capital
                        tempk=k[l]; // the value of capital
                        templ=pow(tempp/eta,eta)*pow(alpha/(1+taxbar),eta-1)*pow((1+phi*tempr)/(tempa*a*pow(tempk,beta)),1/(omega*(eta-1)+alpha)); // the value of labor
                        tempv=pow(tempp/eta,1-alpha/omega)*pow(alpha/(1+taxbar),-alpha/omega)*pow((1+phi*tempr)/(tempa*a*pow(tempk,beta)),omega/(omega*(eta-1)+alpha)); // the value of imported goods
                        tempy=tempa*A*pow(tempk,beta)*pow(templ,alpha)*pow(tempv,eta); // the value of output
                        tempi=(alpha*tempa*A*pow(tempk,beta)*pow(templ,alpha-1)*pow(tempv,eta))/(1+taxbar)-tempp*tempv*(1+phi*tempr); // the value of investment

                        /* conventional maximization step of the Bellman equation (done each "FULL" of "ALL" VF iterations) */
                        kkn[l][j][i]=1; // get ready to start keeping track of where I am in that matrix
                        bbn[l][j][i]=1; // get ready to start keeping track of where I am in that matrix
                        VA[l][j][i]=-999999999; // set the value function to equal to negative infinity
                        for(n=0;n<nb;n=n+1){ // next period's bond
                            for(m=0;m<nk;m=m+1){ // next period's capital
                                tempq=(k[m]-k[l])/k[l]; // the value of the price of equity
                                /* check first if at candidate kkn and bbn the constraint binds, if it does, recalculate L, V, and Y for binding case */
                                if(binding=1){ // if the collateral constraint is binding
                                    tempcolat=-(1-kappa)*(1+a*tempq)*k[m]; // the collateral constraint
                                    if(tempqb*bond[n]-phi*(1+tempr)*(tempp*tempv+(1+taxbar)*pow(templ,omega))<=tempcolat){ // if the collateral constraint is binding
                                        templb=pow((tempqb*bond[n]-tempcolat)/(phi*(1+tempr)*(1+taxbar)*(1+eta/alpha)),1/omega); // labor with collateral constraint binding
                                        tempvb=(tempqb*bond[n]-tempcolat)/(phi*(1+tempr)*tempp*(1+alpha/eta)); // foreign inputs under collateral constraint
                                        if(templb>=0 & tempvb>=0){
                                            tempyb=tempa*A*pow(tempk,beta)*pow(templb,alpha)*pow(tempvb,eta); // output under collateral constraint
                                            p1nval=(tempyb-tempp*tempvb-(k[m]-k[l])*(1+a*tempq/2)-phi*tempr*(tempp*tempvb+(1+taxbar)*pow(templb,omega))-k[l]*delta+bond[j]-tempqb*bond[n])/(1+taxbar)-pow(templb,omega)/omega; // the value in the utility function with the collateral constraint binding
                                        }
                                    } // when the constraint binds, what the value in the utility function is
                                } // if the collateral constraint is binding, use the values of variables under the binding constraint
                                else{
                                    p1nval=(tempy-tempp*tempv-(k[m]-k[l])*(1+a*tempq/2)-phi*tempr*(tempp*tempv+(1+taxbar)*pow(templ,omega))-tempk*delta+bond[j]-tempqb*bond[n])/(1+taxbar)-pow(templ,omega)/omega; // the value in the utility function when the collateral constraint is not binding
                                } // when the constraint doesn't bind, what the value in the utility function is

                                if(p1nval > 0){ // this is when the situation is normal; otherwise we don't consider this situation
                                    expv=0;
                                    for(o=0;o<ne;o=o+1){
                                        if(V[m][n][o]>=-999999999){
                                            expv=pi[o][i]*V[m][n][o]+expv; // p are the transition probabilities, V are the value functions
                                        }
                                    } // get the expected value
                                    p1n[m][n]=pow(p1nval,1-sigma)/(1-sigma)+pow(1+p1nval,-gamma)*expv; //
                                    if(p1n[m][n]>=VA[l][j][i]){
                                        kkn[l][j][i]=m; // keep track of the largest value; which value of capital achieves this
                                        bbn[l][j][i]=n; // keep track of the largest value; which value of bond achieves this
                                        VA[l][j][i]=p1n[m][n];
                                    } // if the value function has a larger value for these values of next period's k and b, then the value function takes this larger value; otherwise, the value function stays the same
                                }
                            } // loop for next period's capital
                        } // loop for next period's bond
                        if(skip==1){ // if V is similar to VA, change to a new k
                            continue;
                        }
                        if(con5<full){ // for the first two of 26 iterations, change to a new k
                            continue;
                        }
//                VA[l][j][i]=(tempy-tempp*tempv-(k[m]-k[l])*(1+a/2*(k[m]-k[l])/k[l])-phi*tempr*(tempp*tempv+(1+taxbar)*pow(templ,omega))-k[l]*delta+bond[j]-bond[n]/(tempr+1))/(1+taxbar)-pow(templ,omega)/omega;
                        if(VA[l][j][i]<0){
                            VA[l][j][i]=-999999999;
                            continue; // pick a new k
                        }
                        expv=0;
                        for(o=0;o<ne;o=o+1){
                            if(V[m][n][o]<-999999999){
                                VA[l][j][i]=-999999999;
                                continue;
                            }
                            expv=pi[o][i]*V[m][n][o]+expv;
                        }
                        VA[l][j][i]=pow(VA[l][j][i],1-sigma)/(1-sigma)+expv/pow(1+VA[l][j][i],gamma);
                    } // loop for this period's capital
                } // loop for this period's bonds
            } // loop for this period's shocks
lifeisgood (51)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
 
/* Testing convergence of decision rules */
    for(i=0;i<nk;i=i+1){
        for(j=0;j<nb;j=j+1){
            for(l=0;l<ne;l=l+1){
                jjkk[i][j][l]=abs(kkn[i][j][l]-kzt[i][j][l])+abs(kkn[i][j][l]-kz2[i][j][l]); // the changes from last period to this period, and from two periods ago to this period
                jjbb[i][j][l]=abs(bbn[i][j][l]-bzt[i][j][l])+abs(bbn[i][j][l]-bz2[i][j][l]); // when we move in the grid of k by b by shock, see if the b' that maximizes VF changes
            }
        }
    }
    for(i=0;i<nk;i=i+1){
        for(j=0;j<nb;j=j+1){
            jjk=jjk+ne-std::count(jjkk[i][j],jjkk[i][j]+8,0); // see whether the k' that maximizes the VF has changed from last period and/or two periods ago; this variable should take the values of 0, 1, or 2
            jjb=jjb+ne-std::count(jjbb[i][j],jjbb[i][j]+8,0); // if the b' is the same in this period, next period and the period after the next, then jjb=0; if it is different from its value of next period and two periods later, then jjb=2
        }
    }
    double rules[nk][nb][ne];
    for(i=0;i<ne;i=i+1){
        for(j=0;j<nb;j=j+1){
            for(l=0;l<nk;l=l+1){
                rules[l][j][i]=abs(V[l][j][i]-VA[l][j][i]);
            }
        }
    }
    bool mask[nk][nb][ne];
    for(i=0;i<ne;i=i+1){
        for(j=0;j<nb;j=j+1){
            for(l=0;l<nk;l=l+1){
                mask[l][j][i]=(rules[l][j][i]<800000);
                if(rules[l][j][i]>=mask[l][j][i]){
                    test2_a=rules[l][j][i];
                }
                else{
                    test2_a=mask[l][j][i];
                }
                if(test2_a>test2){
                    test2=test2_a;
                }
            }
        }
    }
    if(jjk+jjb>=1|skip==0){ // if k' and b' are different than those from the the two previous time periods or if V doesn't equal to VA
        if(con<maxit){ // if we are not at the last iteration
            for(i=0;i<nk;i=i+1){
                for(j=0;j<nb;j=j+1){
                    for(l=0;l<ne;l=l+1){
                        V[i][j][l]=VA[i][j][l]; // update the value of V to prepare for the next iteration
                        VA[i][j][l]=0; // update the value of VA to prepare for the next iteration
                    }
                }
            }
            jjk=nk*nb*ne;
            jjb=jjk;
            if(con5==all){ // if we have reached 26 iterations
                con5=0; // reset
            }
        }
    }
            else if(jjk+jjb==0&skip==1){
                leave=true;
                break; // leave while loop
            }
        } // while loop for whenever con<2000
        if(leave)break; // leave for loop
    } // for loop for every 26 iterations
}
bheadmaster (51)
A friendly suggestion - try to localize the problem and post the smallest possible self-contained code that recreates the error. That way the probability that someone will look over your code and help you with the issue grows exponentially :)
Just the sight of that 572 lines blob of code (which obviously can't even fit a single post!) chilled me to the bone.

Also, localization is the first step of debugging. You might even figure out the error yourself once you break your code down to parts and test them one by one.
Last edited on
jlb (4973)
Be careful, you seem to have quite a few out of bounds array accesses in your program.
1
2
3
4
5
6
7
8
9
10
11
12
 
...
    double pip[2][2]; // transition probabilities between the two states of P
    pip[1][1]=(1-rhop)*pp + rhop;
    pip[2][1]=1-pip[1][1];
    pip[2][2]=(1-rhop)*pp + rhop;
    pip[1][2]=1-pip[2][2];
...
    double piar[4][4]; // transition probabilities for TFP and R combined
...
    double pi[8][8];
...
    pi[8][4]=pip[2][2]*piar[4][2];

You're accessing your arrays out of bounds. Remember that arrays start at zero and end at size - 1, you're accessing array[size] and in some cases you appear to be starting at 1. Remember that with an array of size 2 the valid elements are 0 and 1, using 2 would access the array out of bounds.



lifeisgood (51)
Thank you for the suggestion! My friend has probably also mentioned something like this. Localization means making little .cpp files that contain only a portion of the code? Like I pick out a portion of the code, make a .cpp file of it, compile, and run it?
bheadmaster (51)
Okay, I admit it, I glanced over your code a little bit, and found a fundamental error:

1
2
3
4
5
 
double pi[8][8];
    pi[1][1]=pip[1][1]*piar[1][1];
    pi[2][1]=pip[1][1]*piar[2][1];
    // [...]
    pi[8][4]=pip[2][2]*piar[4][2]; // pi[8] is the 9-th element!!!


Arrays in C++ are 0-indexed, which means that if you create an array of length 8, like you did up there, you can only access indices 0-7. Trying to write to the index 8 causes undefined behaviour.

Try to fix that and see if there any segmentation errors happen then.


EDIT:
Yes, localization means extracting only the part of the code that causes errors.
For example, if you have a program with a dozen of functions called from your main() function, and you figure out that the error happens after the 5-th function is called, and before the 6-th function is called, it's most likely that the 5-th or 6-th function causes error. That makes looking for mistakes a hella lot easier.
Last edited on
jlb (4973)
Also if you ran the program with your debugger the debugger should be able to tell you exactly where it detected the problem. Then you can backtrace to find the actual location of the problem.

And note it's not just pi[][] that is the problem, you're accessing most of your arrays out of bounds in multiple places, multiple times.
lifeisgood (51)
@jlb: I fixed the out of bounds accesses and recompiled, but the debugger is still telling me the same mistake: program received signal SIGSEGV, Segmentation fault. I am now going to search how to have the debugger tell me exactly where the problem is. Right now I am using a graphical editor called Code::Blocks on a Windows computer with GCC compiler.
Last edited on
lifeisgood (51)
@bheadmaster: I fixed the multiple places where I tried to access out of bound, but still get the same error message. How do I figure out where the problem occurs? I am going to Google for this information as well.
jonnin (11329)
do one thing at a time until it crashes by commenting out code, or step through the code line by line in a debugger.

there is a time and a place to unroll a loop, sure, but if you break the index you have a lot of code to muddle through....


if you don't have a good debugger handy just add print statements and comment out chunks until you isolate the issue.
Last edited on
jlb (4973)
Right now I am using a graphical editor called Code::Blocks on a Windows computer with GCC compiler.

Then all you should need to do is run the program with the debugger, F8, or from the menu select Debug then start, or press the red arrow on the debug menu.

I fixed the out of bounds accesses and recompiled

Did you fix them all? There are quite a few out of bounds array access in your program. You need to look at the sizes of your arrays and where your accessing your arrays. You don't seem to understand that arrays start at zero and stop at size - 1.
Last edited on
lifeisgood (51)
@jlb I think so. I was following those guidelines until the lines with multiplication, when I was carried away by the Fortran code that I was imitating. Here's my code at the moment:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
 
/*
* File: main.cpp
* Author: Lia Yin
*
* Created on June 5, 2017
*/

#include <iostream>
#include <vector>
#include <algorithm> // std::count
#include <functional>
#include <tuple>
#include <cmath>

/* Includes, system */
#include <iomanip>
#include <fstream>

// define the variables and parameters
#define nk 60
#define nb 80
#define ne 8
#define na 2
#define nr 2
#define np 2

int skip;

/*
 * used in the section for testing the convergence of decision rules
 */
bool IsNonzero(int i){
    return(i!=0);
} // if the number is not equal to 0, IsNonzero returns a 1

int main(){
    double test2; // criterion for the difference between V and VA
    double test2_a;

    int i; // this one and the ones below are counters
    int j;
    int l;
    int m;
    int n;
    int o;

    // these variables are for the VFs
    double tempa; // the value of TFP for a specific combination of shocks
    double tempr; // the value of world interest rate
    double tempqb; // the value of bond price
    double tempp; // the value of int'l goods price
    double tempk; // the value of capital

    double templ; // the value of labor
    double tempv; // the value of imported goods
    double tempy; // the value of output

    double tempi; // the value of investment

    double tempq; // the value of the price of equity
    double tempcolat; // the collateral constraint

    double templb; // labor with collateral constraint binding
    double tempvb; // foreign inputs under collateral constraint
    double tempyb; // output under collateral constraint

    double p1nval; // value in the parentheses of the utility function
    double p1n[nk][nb]={-999999999}; // for each pair of k' and b', what is the value function equal to

    double expv; // expected value

    int kkn[nk][nb][ne]; // counter for next period's capital
    int bbn[nk][nb][ne]; // counter for next period's bond
    int kzt[nk][nb][ne]; // counter for last period's capital
    int kz2[nk][nb][ne]; // counter for capital from two periods ago
    int bzt[nk][nb][ne]; // counter for last period's bond
    int bz2[nk][nb][ne]; // counter for bond from two periods ago

    int con5; // counter for breaking the iterations down into smaller pieces: in total there are 2000 iterations; break them down into pieces of 26 iterations each
    int con; // counter for total number of iterations; the number is 2000

    int jjkk[nk][nb][ne]; // tally up to see if this group of iterations has the same k' value as last group of iterations and the group of iterations from two times ago
    int jjbb[nk][nb][ne];

    int jjk=0;
    int jjb=0;

    bool binding; // whether the collateral constraint is binding

    double VA[nk][nb][ne]; // updated value after each value function iteration

    bool leave=false; // break out of the loop if necessary


/*
 * INITIAL VALUES FOR PARAMETERS
 */
    const double maxit=2000; // max VF iterations with first price guess
    const double full=2; // VF iterations maximizing over k,b grids out of a total of VF iterations defined in "ALL"
    const double all=26; // "all-full" VF iterations do not maximize, just iterate using VF & dec rules of previous iteration
    double switching=0.00000005; // sets value of "VFDIF" at which VF ITERATIONS SWITCH TO MAXIMIZATION AT ALL TIMES

/* Mexican calibration values */
    const double sigma=2; // CRRA coefficient
    const double rbar=0.08571; // world real interest rate
    const double omega=1.8461; // labor elasticity in the utility function
    const double alpha=0.5924043; // share of capital in the production function
    const double beta=0.3053667; // share of labor in the production function
    const double eta=1-alpha-beta; // share of foreign goods in the production function
    const double phi=0.2579; // fraction of inputs of V and m financed by working capital from firm-level data
    const double A=6.982; // productivity coefficient
    const double kappa=1e-150; // ceiling on the leverage ratio; 1*10^(-149.5701605)
    const double a=2.75; // capital adjustment cost coefficient
    const double delta=0.08800443; // capital depreciation rate
    const double gamma=0.01660445; // semielasticity of the rate of time preference with respect to composite good c-N(L)

/* steady state values */
    const double abar=1; // steady state technology
    const double taxbar=0.169; // steady state tax
    const double pbar=1.028; // steady state price

    const double kbar = 799.592; // steady state capital
    const double lbar = 19.050; // steady state labor
    const double vbar = 45.241; // steady state foreign goods
    const double cbar=265.43703344; // steady state consumption

    double wel=cbar-pow(lbar,omega)/omega; // one component of the steady state utility function
    const double rate=exp(gamma*log(1+wel)); // steady state time preference
    const double wel1=pow(wel,1-sigma)/(1-sigma); // steady state period utility function
    wel=wel1*(1/(1-1/rate)); // steady state total expected utility function

    double V[nk][nb][ne]={wel}; // initial guess of the value function

/* specification of markov process of exogenous shocks */
    double B[nb]; // bond grid with nb elements
    double K[nk]; // capital grid with nk elements
    double E[ne]; // shock grid with ne elements
    double P[ne][ne]={0}; // matrix of transition probabilities, initial values of which are all 0.0
    double Ashock[2]; // vector of TFP shocks
    Ashock[0] = log(1-0.0134);Ashock[1] = log(1+0.0134);
    double Rshock[2]; // vector of int'l rate shocks
    Rshock[0] = log(1-0.0196);Rshock[1] = log(1+0.0196);
    double Pshock[2]; // vector of int'l goods price shocks
    Pshock[0] = log(1-0.0335);Pshock[1] = log(1+0.0335);


/*
 * CALCULATE THE STATE TRANSITION MATRIX
 */

/* A. shocks in int'l goods prices are independent of shocks in TFP or shocks in int'l interest rate, so let's take care of this part first */
    double rhop = 0.737; // autocorrelation coefficient for int'l goods prices
    double pp = 0.5; // if we just consider P, what's the probability its shock is of a certain form
    double pip[2][2]; // transition probabilities between the two states of P
    pip[0][0]=(1-rhop)*pp + rhop;
    pip[1][0]=1-pip[0][0];
    pip[1][1]=(1-rhop)*pp + rhop;
    pip[0][1]=1-pip[1][1];

/* B. joint Markov process for TFP and R shocks */
    double rhoa = 0.537; // autocorrelation coefficient for TFP
    double rhor = 0.572; // autocorrelation coefficient for int'l rate
    double rhoar = (rhoa+rhor)/2; // the procedure requires that the AR(1) coefficients of the shocks that are correlated with each other be the same
    double corrar = -0.67; // correlation coefficient between TFP and int'l rate

    double par[4]; // the probabilities of the four different combinations of TFP and R
    par[0] = (1+corrar)/4; // the probability that both shocks are positive
    par[1]=0.5-par[0];
    par[2]=par[1];
    par[3]=par[2];

    double piar[4][4]; // transition probabilities for TFP and R combined
    int b; // a boolean
    for(i=0;i<4;i=i+1) {
        for(j=0;j<4;j=j+1) {
            if(i=j) {
                b=0;
            }
            else{
                b=1;
            }
            piar[i][j]=(1-rhoar)*par[i]+rhoar*b;
        }
    }

lifeisgood (51)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
 
/* C. Combine all three things, TFP, R, and P, make the large transition probability matrix */
    double pi[8][8];
    pi[0][0]=pip[0][0]*piar[0][0];
    pi[1][0]=pip[0][0]*piar[1][0];
    pi[2][0]=pip[1][0]*piar[0][0];
    pi[3][0]=pip[1][0]*piar[1][0];
    pi[4][0]=pip[0][0]*piar[2][0];
    pi[5][0]=pip[0][0]*piar[3][0];
    pi[6][0]=pip[1][0]*piar[2][0];
    pi[7][0]=pip[1][0]*piar[3][0];
    pi[0][1]=pip[0][0]*piar[0][1];
    pi[1][1]=pip[0][0]*piar[1][1];
    pi[2][1]=pip[1][0]*piar[0][1];
    pi[3][1]=pip[1][0]*piar[1][1];
    pi[4][1]=pip[0][0]*piar[2][1];
    pi[5][1]=pip[0][0]*piar[3][1];
    pi[6][1]=pip[1][0]*piar[2][1];
    pi[7][1]=pip[1][0]*piar[3][1];
    pi[0][2]=pip[0][1]*piar[0][0];
    pi[1][2]=pip[0][1]*piar[1][0];
    pi[2][2]=pip[1][1]*piar[0][0];
    pi[3][2]=pip[1][1]*piar[1][0];
    pi[4][2]=pip[0][1]*piar[2][0];
    pi[5][2]=pip[0][1]*piar[3][0];
    pi[6][2]=pip[1][1]*piar[2][0];
    pi[7][2]=pip[1][1]*piar[3][0];
    pi[0][3]=pip[0][1]*piar[0][1];
    pi[1][3]=pip[0][1]*piar[1][1];
    pi[2][3]=pip[1][1]*piar[0][1];
    pi[3][3]=pip[1][1]*piar[1][1];
    pi[4][3]=pip[0][1]*piar[2][1];
    pi[5][3]=pip[0][1]*piar[3][1];
    pi[6][3]=pip[1][1]*piar[2][1];
    pi[7][3]=pip[1][1]*piar[3][1];

    pi[0][4]=pip[0][0]*piar[0][2];
    pi[1][4]=pip[0][0]*piar[1][1];
    pi[2][4]=pip[1][0]*piar[0][2];
    pi[3][4]=pip[1][0]*piar[1][2];
    pi[4][4]=pip[0][0]*piar[2][2];
    pi[5][4]=pip[0][0]*piar[3][2];
    pi[6][4]=pip[1][0]*piar[2][2];
    pi[7][4]=pip[1][0]*piar[3][2];
    pi[0][5]=pip[0][0]*piar[0][3];
    pi[1][5]=pip[0][0]*piar[1][3];
    pi[2][5]=pip[1][0]*piar[0][3];
    pi[3][5]=pip[1][0]*piar[1][3];
    pi[4][5]=pip[0][0]*piar[2][3];
    pi[5][5]=pip[0][0]*piar[3][3];
    pi[6][5]=pip[1][0]*piar[2][3];
    pi[7][5]=pip[1][0]*piar[3][3];
    pi[0][6]=pip[0][1]*piar[0][2];
    pi[1][6]=pip[0][1]*piar[1][2];
    pi[2][6]=pip[1][1]*piar[0][2];
    pi[3][6]=pip[1][1]*piar[1][2];
    pi[4][6]=pip[0][1]*piar[2][2];
    pi[5][6]=pip[0][1]*piar[3][2];
    pi[6][6]=pip[1][1]*piar[2][2];
    pi[7][6]=pip[1][1]*piar[3][2];
    pi[0][7]=pip[0][1]*piar[0][3];
    pi[1][7]=pip[0][1]*piar[1][3];
    pi[2][7]=pip[1][1]*piar[0][3];
    pi[3][7]=pip[1][1]*piar[1][3];
    pi[4][7]=pip[0][1]*piar[2][3];
    pi[6][7]=pip[0][1]*piar[3][3];
    pi[6][7]=pip[1][1]*piar[2][3];
    pi[7][7]=pip[1][1]*piar[3][3];

/* joint limiting distribution of R,P,TFP shocks */
    double temppi[8][8];
    for(i=0;i<8;i=i+1){
        for(j=0;j<8;j=j+1){
            temppi[i][j]=pi[i][j];
        }
    } // initialize temppi=pi

    double temppi_a[8][8];
    for(i=0;i<300;i=i+1){
        for(j=0;j<8;j=j+1){
            for(m=0;m<8;m=m+1){
                for(l=0;l<8;l=l+1){
                    temppi_a[j][m]=temppi[j][l]*pi[l][m]+temppi_a[j][m];
                    temppi[j][m]=temppi_a[j][m];
                }
            }
        } // finally finished matrix multiplication
        for(j=0;j<8;j=j+1){
            for(l=0;l<8;l=l+1){
                temppi[j][l]=temppi_a[j][l];
            }
        } // change the values of temppi for the next round of multiplication
    } // transition matrix powered 300 times
    double prob[ne];
    for(i=0;i<ne;i=i+1){
        prob[i]=temppi[i][1]; // becomes a row vector
    } // joint limiting probability of shocks
    m=0;
    double sa[8];
    double sr[8];
    double qb[8];
    double sp[8];
    for(i=0;i<nr;i=i+1){
        for(j=0;j<np;j=j+1){
            for(l=0;l<na;l=l+1){
                m=m+1;
                sa[m]=exp(Ashock[l])*abar; // state for TFP
                sr[m]=exp(Rshock[i])*(1+rbar); // state for the interest rate
                qb[m]=1/sr[m]; // state for bond prices
                sp[m]=exp(Pshock[j])*pbar; // state for price of int'l goods
            }
        }
    }

/*
 * CREATE GRIDS
 */

/* A. domestic capital */
    double kss=752.270;
    double kgridgap=3.2;
    double klow=kss-kgridgap*(nk/2);
    double k[nk];
    for(i=0;i<nk;i=i+1){
        k[i]=klow+i*kgridgap;
    }

/* B. foreign bonds */
    double bss=-127.184;
    double bgridgap=2*kgridgap;
    double blow=bss-bgridgap*(nb/2);
    double bond[nb];
    for(i=0;i<nb;i=i+1){
        bond[i]=blow+i*bgridgap;
    }


/*
 * CALCULATE THE VALUE OF CHOOSING CAPITAL K'(M) AND ASSET HOLDINGS B'(N) WHEN THE INITIAL STATE IS
 * GIVEN BY (L,J,I)
 */

    skip=0;
    test2=switching;
lifeisgood (51)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
 
/* A. calculate labor, int'l inputs allocations for given K(I), S(IS) assuming credit constraint does
 not bind; note that they are independent of K' and B' */

    for(i=0;i<ne;i=i+1){
        for(j=0;j<nb;j=j+1){
            for(l=0;l<nk;l=l+1){
                kzt[i][j][l]=-999999999;
                kz2[i][j][l]=-999999999;
            }
        }
    } // initialize the values of kzt and kz2 to make sure that we don't stop at the first iteration

    for(con5=0;con5<all;con5=con5+1){ // every 26 iterations
        con=con+1; // add one every time it loops
        while(con<maxit){ // whenever con<2000
            if(test2<switching){
                skip=1; // when V is similar to VA, restart maximization
            }
            for(i=0;i<ne;i=i+1){ // this period's shocks
                tempa=sa[i]; // the value of TFP for this specific combination of shocks
                tempr=sr[i]-1; // the value of world interest rate
                tempqb=qb[i]; // the value of bond price
                tempp=sp[i]; // the value of int'l goods price
                for(j=0;j<nb;j=j+1){ // this period's bond
                    for(l=0;l<nk;l=l+1){ // this period's capital
                        tempk=k[l]; // the value of capital
                        templ=pow(tempp/eta,eta)*pow(alpha/(1+taxbar),eta-1)*pow((1+phi*tempr)/(tempa*a*pow(tempk,beta)),1/(omega*(eta-1)+alpha)); // the value of labor
                        tempv=pow(tempp/eta,1-alpha/omega)*pow(alpha/(1+taxbar),-alpha/omega)*pow((1+phi*tempr)/(tempa*a*pow(tempk,beta)),omega/(omega*(eta-1)+alpha)); // the value of imported goods
                        tempy=tempa*A*pow(tempk,beta)*pow(templ,alpha)*pow(tempv,eta); // the value of output
                        tempi=(alpha*tempa*A*pow(tempk,beta)*pow(templ,alpha-1)*pow(tempv,eta))/(1+taxbar)-tempp*tempv*(1+phi*tempr); // the value of investment

                        /* conventional maximization step of the Bellman equation (done each "FULL" of "ALL" VF iterations) */
                        kkn[l][j][i]=1; // get ready to start keeping track of where I am in that matrix
                        bbn[l][j][i]=1; // get ready to start keeping track of where I am in that matrix
                        VA[l][j][i]=-999999999; // set the value function to equal to negative infinity
                        for(n=0;n<nb;n=n+1){ // next period's bond
                            for(m=0;m<nk;m=m+1){ // next period's capital
                                tempq=(k[m]-k[l])/k[l]; // the value of the price of equity
                                /* check first if at candidate kkn and bbn the constraint binds, if it does, recalculate L, V, and Y for binding case */
                                if(binding=1){ // if the collateral constraint is binding
                                    tempcolat=-(1-kappa)*(1+a*tempq)*k[m]; // the collateral constraint
                                    if(tempqb*bond[n]-phi*(1+tempr)*(tempp*tempv+(1+taxbar)*pow(templ,omega))<=tempcolat){ // if the collateral constraint is binding
                                        templb=pow((tempqb*bond[n]-tempcolat)/(phi*(1+tempr)*(1+taxbar)*(1+eta/alpha)),1/omega); // labor with collateral constraint binding
                                        tempvb=(tempqb*bond[n]-tempcolat)/(phi*(1+tempr)*tempp*(1+alpha/eta)); // foreign inputs under collateral constraint
                                        if(templb>=0 & tempvb>=0){
                                            tempyb=tempa*A*pow(tempk,beta)*pow(templb,alpha)*pow(tempvb,eta); // output under collateral constraint
                                            p1nval=(tempyb-tempp*tempvb-(k[m]-k[l])*(1+a*tempq/2)-phi*tempr*(tempp*tempvb+(1+taxbar)*pow(templb,omega))-k[l]*delta+bond[j]-tempqb*bond[n])/(1+taxbar)-pow(templb,omega)/omega; // the value in the utility function with the collateral constraint binding
                                        }
                                    } // when the constraint binds, what the value in the utility function is
                                } // if the collateral constraint is binding, use the values of variables under the binding constraint
                                else{
                                    p1nval=(tempy-tempp*tempv-(k[m]-k[l])*(1+a*tempq/2)-phi*tempr*(tempp*tempv+(1+taxbar)*pow(templ,omega))-tempk*delta+bond[j]-tempqb*bond[n])/(1+taxbar)-pow(templ,omega)/omega; // the value in the utility function when the collateral constraint is not binding
                                } // when the constraint doesn't bind, what the value in the utility function is

                                if(p1nval > 0){ // this is when the situation is normal; otherwise we don't consider this situation
                                    expv=0;
                                    for(o=0;o<ne;o=o+1){
                                        if(V[m][n][o]>=-999999999){
                                            expv=pi[o][i]*V[m][n][o]+expv; // p are the transition probabilities, V are the value functions
                                        }
                                    } // get the expected value
                                    p1n[m][n]=pow(p1nval,1-sigma)/(1-sigma)+pow(1+p1nval,-gamma)*expv; //
                                    if(p1n[m][n]>=VA[l][j][i]){
                                        kkn[l][j][i]=m; // keep track of the largest value; which value of capital achieves this
                                        bbn[l][j][i]=n; // keep track of the largest value; which value of bond achieves this
                                        VA[l][j][i]=p1n[m][n];
                                    } // if the value function has a larger value for these values of next period's k and b, then the value function takes this larger value; otherwise, the value function stays the same
                                }
                            } // loop for next period's capital
                        } // loop for next period's bond
                        if(skip==1){ // if V is similar to VA, change to a new k
                            continue;
                        }
                        if(con5<full){ // for the first two of 26 iterations, change to a new k
                            continue;
                        }
//                VA[l][j][i]=(tempy-tempp*tempv-(k[m]-k[l])*(1+a/2*(k[m]-k[l])/k[l])-phi*tempr*(tempp*tempv+(1+taxbar)*pow(templ,omega))-k[l]*delta+bond[j]-bond[n]/(tempr+1))/(1+taxbar)-pow(templ,omega)/omega;
                        if(VA[l][j][i]<0){
                            VA[l][j][i]=-999999999;
                            continue; // pick a new k
                        }
                        expv=0;
                        for(o=0;o<ne;o=o+1){
                            if(V[m][n][o]<-999999999){
                                VA[l][j][i]=-999999999;
                                continue;
                            }
                            expv=pi[o][i]*V[m][n][o]+expv;
                        }
                        VA[l][j][i]=pow(VA[l][j][i],1-sigma)/(1-sigma)+expv/pow(1+VA[l][j][i],gamma);
                    } // loop for this period's capital
                } // loop for this period's bonds
            } // loop for this period's shocks
Last edited on
lifeisgood (51)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
 
/* Testing convergence of decision rules */
    for(i=0;i<nk;i=i+1){
        for(j=0;j<nb;j=j+1){
            for(l=0;l<ne;l=l+1){
                jjkk[i][j][l]=abs(kkn[i][j][l]-kzt[i][j][l])+abs(kkn[i][j][l]-kz2[i][j][l]); // the changes from last period to this period, and from two periods ago to this period
                jjbb[i][j][l]=abs(bbn[i][j][l]-bzt[i][j][l])+abs(bbn[i][j][l]-bz2[i][j][l]); // when we move in the grid of k by b by shock, see if the b' that maximizes VF changes
            }
        }
    }
    for(i=0;i<nk;i=i+1){
        for(j=0;j<nb;j=j+1){
            jjk=jjk+ne-std::count(jjkk[i][j],jjkk[i][j]+8,0); // see whether the k' that maximizes the VF has changed from last period and/or two periods ago; this variable should take the values of 0, 1, or 2
            jjb=jjb+ne-std::count(jjbb[i][j],jjbb[i][j]+8,0); // if the b' is the same in this period, next period and the period after the next, then jjb=0; if it is different from its value of next period and two periods later, then jjb=2
        }
    }
    double rules[nk][nb][ne];
    for(i=0;i<ne;i=i+1){
        for(j=0;j<nb;j=j+1){
            for(l=0;l<nk;l=l+1){
                rules[l][j][i]=abs(V[l][j][i]-VA[l][j][i]);
            }
        }
    }
    bool mask[nk][nb][ne];
    for(i=0;i<ne;i=i+1){
        for(j=0;j<nb;j=j+1){
            for(l=0;l<nk;l=l+1){
                mask[l][j][i]=(rules[l][j][i]<800000);
                if(rules[l][j][i]>=mask[l][j][i]){
                    test2_a=rules[l][j][i];
                }
                else{
                    test2_a=mask[l][j][i];
                }
                if(test2_a>test2){
                    test2=test2_a;
                }
            }
        }
    }
    if(jjk+jjb>=1|skip==0){ // if k' and b' are different than those from the the two previous time periods or if V doesn't equal to VA
        if(con<maxit){ // if we are not at the last iteration
            for(i=0;i<nk;i=i+1){
                for(j=0;j<nb;j=j+1){
                    for(l=0;l<ne;l=l+1){
                        V[i][j][l]=VA[i][j][l]; // update the value of V to prepare for the next iteration
                        VA[i][j][l]=0; // update the value of VA to prepare for the next iteration
                    }
                }
            }
            jjk=nk*nb*ne;
            jjb=jjk;
            if(con5==all){ // if we have reached 26 iterations
                con5=0; // reset
            }
        }
    }
            else if(jjk+jjb==0&skip==1){
                leave=true;
                break; // leave while loop
            }
        } // while loop for whenever con<2000
        if(leave)break; // leave for loop
    } // for loop for every 26 iterations
}
Last edited on
lifeisgood (51)
@jonnin Those are really good suggestions. Right now I am not able to step the code in a debugger because it crashes as soon as I start the debugger, but I am looking into the print statement possibility.
MikeyBoy (5631)
If the program crashes when you run it in the debugger, then you can use the debugger to look at the call stack at the time of the crash. That should tell you where in the code the crash is happening.
jlb (4973)
You really really need to learn to use functions.

And you need to start using more meaningful variable names (names that actually describe their purpose). And you really should avoid using variables like 'l' (lower case L), 'O' (upper case o) as these variables can be hard to distingush between 1 (one) and 0 (Zero). You also seem to have forgot that the operator= is the assignment operator, your using it in places where you should be using the comparison operator==.

If it is crashing as soon as you start the debugger the debugger should point you to the line where the problem is detected.

Did you try to set a breakpoint early in the program and then single step through the program?

If you can't figure out the debugger you should try putting output statements at various parts of the program so you can start to tell how far the program is running before it crashes. Since your only function is so large I suggest you put a output statement at the first line, then probably about every 30 lines there after. All the output needs to be at this point is something like "I'm at line 1", with the number reflecting the current line and be sure to flush the output after every output.

Pages: 123