Оценивание линейных регрессионных моделей в условиях гетероскедастичности возмущений., страница 3

#include<math.h>

#include <stdlib.h>

#include<iostream.h>

#include <fstream.h>

#include<io.h>

const k=2; //number factors

const m=5; //number model

const n=36;//number experiments  6x6

/********* teta for model ********/

float teta[]={ 1.0, 0.01, 1.0, 0.1, 1.0};

/********* drop factors **********/

float r[]={-1.0,-0.6,-0.2,0.2,0.6,1.0};

/* function from factors of exit */

float f1(float x1,float x2)

{  return 1/x1;         }

float f2(float x1,float x2)

{  return x1;           }

float f3(float x1,float x2)

{  return x2*x2;        }

float f4(float x1,float x2)

{  return x2;           }

float f5(float x1,float x2)

{  return 1;            }

/********* element eta  **********/

float el_eta(float x1,float x2,float *t)

{ return(t[0]*f1(x1,x2)+t[1]*f2(x1,x2)+t[2]*f3(x1,x2)+t[3]*f4(x1,x2)+t[4]*f5(x1,x2));}

/********* dispersion ************/

float sigma2(float u[n])

{ float ro,sigma2_T,usr=0;

int i;

for (i=0; i<n; i++) usr=usr+u[i];

usr=usr/n;

sigma2_T=0;

for(i=0;i<n;i++) sigma2_T=sigma2_T+(u[i]-usr)*(u[i]-usr);

sigma2_T=sigma2_T/(n-1);          // D[u]

ro=0.05;                          // 5%

return(ro*ro*sigma2_T);

}

/******** mod from a ****/

float F(float a)

{ while (a>0) a=a-1;

return(a+1);

}

/******** randomize ****/

void psld_ksi(float*u)

{ float ksi,pi=M_PI;

ksi=F(11*0.5+pi);

u[0]=2*ksi-1;

for(int i=1;i<n;i++)

{ ksi=F(11*ksi+pi);

u[i]=2*ksi-1;

}

return;

}

/*************************/

void matrix(float a[n][m])

{int i,j,k;

k=0;

for(i=0;i<6;i++)

for(j=0;j<6;j++)

{ a[k][0]=f1(r[i],r[j]);

a[k][1]=f2(r[i],r[j]);

a[k][2]=f3(r[i],r[j]);

a[k][3]=f4(r[i],r[j]);

a[k][4]=f5(r[i],r[j]);

k++;

}

}

//_________________________________//

void lab1(float u[n],float eps[n],float y[n])

{int i,k,j;

float sg;

psld_ksi(u);

sg=sqrt(sigma2(u));

cout<<"sigma true = "<<sg*sg<<"\n";

for(i=0;i<n;i++)

{ float s;

s=u[i]*u[i]+u[i+1]*u[i+1];

if(s<1) eps[i]=sg*u[i]*sqrt(-2*log(s)/s);

}

k=0;

for(i=0;i<6;i++)

for(j=0;j<6;j++)

{ y[k]=el_eta(r[i],r[j],teta)+eps[k];

k++;

}

ofstream out1("y.dat");

for(i=0;i<n;i++) out1<<y[i]<<"\n";

out1.close();

}

////////////////////////////////////////////////////

int min_el(int num,float y[n])

{float a;

int i,j=num;

a=y[num];

for(i=num+1;i<n;i++)

if(a>y[i]) { a=y[i]; j=i; }

return j;

}

/***********************************/

void getera(float y[n],float u[n])

{int i,j;

float mel;

for(i=0;i<n;i++) u[i]=fabs(y[i]);

for(i=0;i<n-1;i++)

{ j=min_el(i,u);

mel=u[i];   u[i]=u[j];    u[j]=mel;

}

for(i=0;i<n;i++) u[i]=0.05*0.05*u[i];  //5%

ofstream out2("V.dat");

for(i=0;i<n;i++) out2<<u[i]<<"\n";

out2.close();

}

/***** scalar multax matrixs *********/

void TXvX(float a[n][m],float u[n],float matr[m][m])

{  for(int i=0;i<m;i++)

for(int j=0;j<m;j++)

{ matr[i][j]=0;

for (int l=0;l<n;l++) matr[i][j]=matr[i][j]+a[l][i]/u[l]*a[l][j];

}

}

/******** scalar multax matrix on vector */

void TXvy(float a[n][m],float y[n],float u[n],float vect[m])

{ for (int i=0;i<m;i++)

{ vect[i]=0;

for (int j=0;j<n;j++) vect[i]=vect[i]+a[j][i]/u[j]*y[j];

}

}

void scal1(float a[n][m],float t[m],float y1[n])

{ for (int i=0;i<n;i++)

{ y1[i]=0;