Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision

Target

Select target project
  • giovanma/dynare
  • giorgiomas/dynare
  • Vermandel/dynare
  • Dynare/dynare
  • normann/dynare
  • MichelJuillard/dynare
  • wmutschl/dynare
  • FerhatMihoubi/dynare
  • sebastien/dynare
  • lnsongxf/dynare
  • rattoma/dynare
  • CIMERS/dynare
  • FredericKarame/dynare
  • SumuduK/dynare
  • MinjeJeon/dynare
  • camilomrch/dynare
  • DoraK/dynare
  • avtishin/dynare
  • selma/dynare
  • claudio_olguin/dynare
  • jeffjiang07/dynare
  • EthanSystem/dynare
  • stepan-a/dynare
  • wjgatt/dynare
  • JohannesPfeifer/dynare
  • gboehl/dynare
  • chskcau/dynare-doc-fixes
27 results
Select Git revision
Show changes
Show whitespace changes
Inline Side-by-side
Showing
with 0 additions and 2186 deletions
dynare++/tests/sw_euro.mod deleted 100644 → 0
View file @ 219d2bb3
var MC EH EF R_KF QF CF IF YF LF PIEF WF RF R_KH QH CH IH YH LH PIEH WH RH EE_A PIE_BAR EE_B EE_G EE_L EE_I KF KH ONE;
varexo E_A E_B E_G E_L E_I ETA_R E_PIE_BAR ETA_Q ETA_P ETA_W ;
parameters xi_e lambda_w alpha czcap beta phi_i tau sig_c hab ccs cinvs phi_y gamma_w xi_w gamma_p xi_p sig_l r_dpi r_pie r_dy r_y rho rho_a rho_pb rho_b rho_g rho_l rho_i ;
alpha=.30;
beta=0.99;
tau=0.025;
ccs=0.6;
cinvs=.22;
lambda_w = 0.5;
phi_i= 6.771;
sig_c= 1.353;
hab= 0.573;
xi_w= 0.737;
sig_l= 2.400;
xi_p= 0.908;
xi_e= 0.599;
gamma_w= 0.763;
gamma_p= 0.469;
czcap= 0.169;
phi_y= 1.408;
r_pie= 1.684;
r_dpi= 0.14;
rho= 0.961;
r_y= 0.099;
r_dy= 0.159;
rho_a= 0.823;
rho_b= 0.855;
rho_g= 0.949;
rho_l= 0.889;
rho_i= 0.927;
rho_pb= 0.924;
model;
CF = (1/(1+hab))*(CF(1)+hab*CF(-1))-((1-hab)/((1+hab)*sig_c))*(RF-PIEF(1)-EE_B) ;
0 = alpha*R_KF+(1-alpha)*WF -EE_A ;
PIEF = 0*ONE;
IF = (1/(1+beta))* (( IF(-1) + beta*(IF(1)))+(1/phi_i)*QF)+0*ETA_Q+EE_I ;
QF = -(RF-PIEF(1))+(1-beta*(1-tau))*((1+czcap)/czcap)*R_KF(1)+beta*(1-tau)*QF(1) +0*EE_I ;
KF = (1-tau)*KF(-1)+tau*IF(-1) ;
YF = (ccs*CF+cinvs*IF)+EE_G ;
YF = 1*phi_y*( alpha*KF+alpha*(1/czcap)*R_KF+(1-alpha)*LF+EE_A ) ;
WF = (sig_c/(1-hab))*(CF-hab*CF(-1)) + sig_l*LF - EE_L ;
LF = R_KF*((1+czcap)/czcap)-WF+KF ;
EF = EF(-1)+EF(1)-EF+(LF-EF)*((1-xi_e)*(1-xi_e*beta)/(xi_e));
CH = (hab/(1+hab))*CH(-1)+(1/(1+hab))*CH(1)-((1-hab)/((1+hab)*sig_c))*(RH-PIEH(1)-EE_B) ;
IH = (1/(1+beta))* (( IH(-1) + beta*(IH(1)))+(1/phi_i)*QH )+1*ETA_Q+1*EE_I ;
QH = -(RH-PIEH(1))+(1-beta*(1-tau))*((1+czcap)/czcap)*R_KH(1)+beta*(1-tau)*QH(1) +EE_I*0+0*ETA_Q ;
KH = (1-tau)*KH(-1)+tau*IH(-1) ;
YH = (ccs*CH+cinvs*IH)+ EE_G ;
YH = phi_y*( alpha*KH+alpha*(1/czcap)*R_KH+(1-alpha)*LH ) +phi_y*EE_A ;
PIEH = (1/(1+beta*gamma_p))*
(
(beta)*(PIEH(1)) +(gamma_p)*(PIEH(-1))
+((1-xi_p)*(1-beta*xi_p)/(xi_p))*(MC)
) + ETA_P ;
MC = alpha*R_KH+(1-alpha)*WH -EE_A;
WH = (1/(1+beta))*(beta*WH(+1)+WH(-1))
+(beta/(1+beta))*(PIEH(+1))
-((1+beta*gamma_w)/(1+beta))*(PIEH)
+(gamma_w/(1+beta))*(PIEH(-1))
-(1/(1+beta))*(((1-beta*xi_w)*(1-xi_w))/(((1+(((1+lambda_w)*sig_l)/(lambda_w))))*xi_w))*(WH-sig_l*LH-(sig_c/(1-hab))*(CH-hab*CH(-1))+EE_L)
+ETA_W;
LH = R_KH*((1+czcap)/czcap)-WH+KH ;
RH = r_dpi*(PIEH-PIEH(-1))
+(1-rho)*(r_pie*(PIEH(-1)-PIE_BAR)+r_y*(YH-YF))
+r_dy*(YH-YF-(YH(-1)-YF(-1)))
+rho*(RH(-1)-PIE_BAR)
+PIE_BAR
+ETA_R;
EH = EH(-1)+EH(1)-EH+(LH-EH)*((1-xi_e)*(1-xi_e*beta)/(xi_e));
EE_A = (rho_a)*EE_A(-1) + E_A;
PIE_BAR = rho_pb*PIE_BAR(-1)+ E_PIE_BAR ;
EE_B = rho_b*EE_B(-1) + E_B ;
EE_G = rho_g*EE_G(-1) + E_G ;
EE_L = rho_l*EE_L(-1) + E_L ;
EE_I = rho_i*EE_I(-1) + E_I ;
ONE = 0*ONE(-1) ;
end;
vcov = [0.357604 0 0 0 0 0 0 0 0 0;
0 0.112896 0 0 0 0 0 0 0 0;
0 0 0.105625 0 0 0 0 0 0 0;
0 0 0 12.39040 0 0 0 0 0 0;
0 0 0 0 0.722500 0 0 0 0 0;
0 0 0 0 0 0.656100 0 0 0 0;
0 0 0 0 0 0 0.000289 0 0 0;
0 0 0 0 0 0 0 0.364816 0 0;
0 0 0 0 0 0 0 0 0.025600 0;
0 0 0 0 0 0 0 0 0 0.083521];
order = 1;
dynare++/tests/swma_pie.dyn deleted 100644 → 0
View file @ 219d2bb3
// this model has sticky wages and adjustment costs in
// investment, consumer goods sector is perfectly competitive, thus MC=1
// with money and transaction costs based on money velocity
// and it has a financial accelerator
// wage is indexed to past consumer price inflation
// LAMBDA Lagrange multiplier on household's budget constraint (divided by price level)
// PIE inflation of CPI
// PIETILDE to what inflation new wage setters index (here PIE(-1) but could be PIEW(-1))
// INT nominal interest rate
// C real consumption
// I real investment
// K real capital
// R real rental rate of capital
// W real wage
// L labour
// Y real output
// PIEW nominal wage inflation
// VW wage front loading term for newly set wages
// BBD, BBE, BBF, BBG terms in nominator and denominator in wage FOC
// G government
// SL process for labor shock
// SC process for consumption shock
// SY process for technology shock
// RM real money balances hold
// Q real price of capital
// Q_M1 lagged Q
// RK nominal return of capital for enterpreneurs
// OMEGABAR threshold value for idiosyncratic shock
// N real net worth of borrowers
// WF lifetime utility
var LAMBDA PIE PIETILDE INT C I K R W L Y PIEW VW BBD BBE BBF BBG G SL SC SY RM
Q Q_M1 RK OMEGABAR N ACAL ACALPRIME BCAL BCALPRIME WF;
varexo E_C E_L E_Y E_GOV E_INT;
parameters dep beta gamma eta biga alpha sigmaw phiw deltaw sg pietar h psi nu osigma mu tc1 tc2 ksi1 ksi2 c_weight rho_g rho_l rho_c rho_y;
dep = 0.025;
beta = 0.99;
gamma = 1;
eta = 2;
alpha = 0.30;
biga = alpha^(-alpha)*(1-alpha)^(alpha-1);
sigmaw = 11;
phiw = 2;
deltaw = 0.75;
sg = 0.18;
pietar = 1.03^0.25;
h = 0.8;
// investment adjustment costs
psi = 12;
// enterpreneur saving rate
nu = 0.94;
// stderr of enterpreneur's idiosyncratic shocks
osigma = 0.5;
// monitoring cost for lender
mu = 0.2;
// consumption transaction costs
tc1 = 0.05;
tc2 = 0.5;
// Taylor rule
ksi1 = 0.106;
ksi2 = 3;
rho_g = 0.90;
rho_l = 0.90;
rho_c = 0.90;
rho_y = 0.90;
// weight of consumption utility
c_weight = 1;
model;
// capital accumulation
K = (1 - dep - psi/2*(I(-1)/K(-1)-dep)^2)*K(-1) + I(-1);
// FOC bonds
LAMBDA = beta*INT*LAMBDA(+1)/PIE(+1);
// FOC consumption (right hand side is equal to LAMBDA*(1+TC+TCPRIME*C/RM))
SC*c_weight*(C-h*C(-1))^(-eta) = LAMBDA*(1+2*tc1*C/RM-2*sqrt(tc1*tc2));
// FOC money (right hand side is equal to 1 - TCPRIME*C*C/RM/RM)
beta*LAMBDA(+1)/LAMBDA/PIE(+1) = 1 - tc1*C*C/RM/RM + tc2;
// FOC investment removed
// FOC capital(+1) removed
// real price of capital
Q = (1-psi*(I/K-dep))^(-1);
// nominal return on capital
RK = PIE*(R + Q*(1 - dep + psi*(I/K-dep)*I/K -psi/2*(I/K-dep)^2))/Q(-1);
// FOC in optimal contract for K(+1)
RK(+1)*(BCAL(+1)*ACALPRIME(+1)/BCALPRIME(+1)-ACAL(+1)) = INT(+1)*ACALPRIME(+1)/BCALPRIME(+1);
// Participation constraint
//RK(+1)*BCAL(+1) = INT(+1)*(1-N(+1)*PIE(+1)/Q/K(+1));
RK*BCAL = INT*(1-N*PIE/Q(-1)/K);
// evolution of net worth (real)
N*PIE*PIE(-1) = nu*(ACAL(-1)+BCAL(-1))*RK(-1)*Q_M1(-1)*K(-1) - nu*INT(-1)*(Q_M1(-1)*K(-1)-N(-1)*PIE);
// marginal cost is 1
1 = biga*(W/SY)^(1-alpha)*R^alpha;
// labor attaining minimal MC
L = (1-alpha)/W*Y;
// capital attaining minimal MC
K = alpha/R*Y;
// FOC for newly set wages
W*VW = sigmaw/(sigmaw-1)*(BBD*VW^(-sigmaw*gamma) + phiw*BBE*VW^(-sigmaw) - phiw*BBF)/BBG;
// definition of BBD
BBD = SL*L^(1+gamma) + deltaw*beta*(PIETILDE(+1)/PIEW(+1))^(-sigmaw*(1+gamma))*BBD(+1);
// definition of BBE
BBE = LAMBDA*L*W + deltaw*beta*(PIETILDE(+1)/PIEW(+1))^(-2*sigmaw)*BBE(+1);
// definition of BBF
BBF = LAMBDA*L*W + deltaw*beta*(PIETILDE(+1)/PIEW(+1))^(-sigmaw)*BBF(+1);
// definition of BBG
BBG = LAMBDA*L + deltaw*beta*(PIETILDE(+1)/PIEW(+1))^(-sigmaw)*PIETILDE(+1)/PIE(+1)*BBG(+1);
// price index
1 = (1-deltaw)*VW^(1-sigmaw) + deltaw*(PIETILDE/PIEW)^(1-sigmaw);
// definition of ACAL
ACAL = 0.5*erfc((log(OMEGABAR) - 0.5*osigma^2)/osigma/sqrt(2.0)) - OMEGABAR/2*erfc((log(OMEGABAR) + 0.5*osigma^2)/osigma/sqrt(2.0));
// definition of BCAL
BCAL = OMEGABAR/2*erfc((log(OMEGABAR) + 0.5*osigma^2)/osigma/sqrt(2.0)) + (1-mu)/2*(1+erf((log(OMEGABAR) - 0.5*osigma^2)/osigma/sqrt(2.0)));
// definition of ACALPRIME
ACALPRIME = -0.5*erfc((log(OMEGABAR) + 0.5*osigma^2)/osigma/sqrt(2.0));
// definition of BCALPRIME
BCALPRIME = -ACALPRIME - mu/osigma/2.506628274631*exp(-((log(OMEGABAR) + 0.5*osigma)^2)/2/osigma/osigma);
// identity for PIEW
PIEW = PIE*W/W(-1);
// welfare identity
WF = SC*c_weight*(C-h*C(-1))^(1-eta)/(1-eta) - SL*L^(1+gamma)/(1+gamma) + beta*WF(+1);
// interest rate rule
INT = INT(-1)^ksi1*((PIE/beta)*(PIE/pietar)^ksi2)^(1-ksi1)*exp(E_INT);
// aggregate constraint
Y = C + I + G + (1-ACAL-BCAL)*RK*Q(-1)*K;
//Y = C + I + G;
// process for government
G/Y = (G(-1)/Y(-1))^rho_g*sg^(1-rho_g)*exp(E_GOV/sg);
// to what do they index (pietar, past inflation, past indexed inflation)
PIETILDE = PIE(-1);
//PIETILDE = pietar;
// exo processes
SL = SL(-1)^rho_l*exp(E_L);
SC = SC(-1)^rho_c*exp(E_C);
SY = SY(-1)^rho_y*exp(E_Y);
// lagged Q
Q_M1 = Q(-1);
end;
initval;
RM = 0.1;
INT = pietar/beta;
PIE = pietar;
PIEW = pietar;
PIETILDE = pietar;
//R = dep/beta;
R = 0.1;
W = (1/biga/(R)^alpha)^(1/(1-alpha));
LAMBDA = ((1-dep*alpha/R-sg)*(1-h)*c_weight/(1-alpha)*W^(1/gamma+1)*((sigmaw-1)/sigmaw)^(1/gamma))^(-1/(1/eta+1/gamma));
L = (W*LAMBDA*(sigmaw-1)/sigmaw)^(1/gamma);
Y = W*L/(1-alpha);
K = alpha/R*Y;
I = dep*K;
G = sg*Y;
VW = 1;
BBD = L^(1+gamma)/(1-deltaw*beta);
BBE = LAMBDA*L*W/(1-deltaw*beta);
BBF = LAMBDA*L*W/(1-deltaw*beta);
BBG = LAMBDA*L/(1-deltaw*beta);
Q = 1;
Q_M1 = Q;
RK = 1/Q*PIE*(R+(1-dep)*Q);
OMEGABAR = 0.5;
ACAL = 0.5*erfc((log(OMEGABAR) - 0.5*osigma^2)/osigma/sqrt(2.0)) - OMEGABAR/2*erfc((log(OMEGABAR) + 0.5*osigma^2)/osigma/sqrt(2.0));
BCAL = OMEGABAR/2*erfc((log(OMEGABAR) + 0.5*osigma^2)/osigma/sqrt(2.0)) + (1-mu)/2*(1+erf((log(OMEGABAR) - 0.5*osigma^2)/osigma/sqrt(2.0)));
ACALPRIME = -0.5*erfc((log(OMEGABAR) + 0.5*osigma^2)/osigma/sqrt(2.0));
BCALPRIME = -ACALPRIME - mu/osigma/2.506628274631*exp(-((log(OMEGABAR) + 0.5*osigma)^2)/2/osigma/osigma);
N = (nu*(ACAL+BCAL)*RK*Q*K-nu*INT*Q*K)/(PIE*PIE-nu*INT*PIE);
C = Y - I - G - (1-ACAL-BCAL)*RK*Q*K;
SL = 1;
SC = 1;
SY = 1;
WF = 1/(1-beta)*(SC*c_weight*((1-h)*C)^(1-eta)/(1-eta) - SL*L^(1+gamma)/(1+gamma));
end;
vcov = [
0.0001 0 0 0 0;
0 0.0001 0 0 0;
0 0 0.0001 0 0;
0 0 0 0.0001 0;
0 0 0 0 0.0001
];
order = 4;
dynare++/tests/test.mod deleted 100644 → 0
View file @ 219d2bb3
var a, b, c, h, k, y;
varexo e,u;
parameters beta, rho, alpha, delta, theta, psi, tau, phi;
alpha = 0.36;
rho = 0.95;
tau = 0.025;
beta = 0.99;
delta = 0.025;
psi = 0;
theta = 2.95;
phi = 0.1;
model;
c*theta*h^(1+psi)=(1-alpha)*y;
k = beta*(((exp(b)*c)/(exp(b(+1))*c(+1)))*(exp(b(+1))*alpha*y(+1)+(1-delta)*k));
y = exp(a)*(k(-1)^alpha)*(h^(1-alpha));
k = exp(b)*(y-c)+(1-delta)*k(-1);
a = rho*a(-1)+tau*b(-1) + e;
b = tau*a(-1)+rho*b(-1) + u;
end;
initval;
y = 1;
c = 0.7;
h = 0.1;
k = 11;
a = 0;
b = 0;
e = 0;
u = 0;
end;
vcov = [ 0.000081 0.000008;0.000008 0.000081];
order = 2;
dynare++/tests/test1.mod deleted 100644 → 0
View file @ 219d2bb3
var y, c, k, a, h, b;
varexo e,u;
parameters beta, rho, alpha, delta, theta, psi, tau, phi;
alpha = 0.36;
rho = 0.95;
tau = 0.025;
beta = 0.99;
delta = 0.025;
psi = 0;
theta = 2.95;
phi = 0.1;
model;
c*theta*h^(1+psi)=(1-alpha)*y;
k = beta*(((exp(b)*c)/(exp(b(+1))*c(+1)))*(exp(b(+1))*alpha*y(+1)+(1-delta)*k));
y = exp(a)*(k(-1)^alpha)*(h^(1-alpha));
k = exp(b)*(y-c)+(1-delta)*k(-1);
a = rho*a(-1)+tau*b(-1) + e;
b = tau*a(-1)+rho*b(-1) + u;
end;
initval;
y = 1;
c = 0.7;
h = 0.1;
k = 11;
a = 0;
b = 0;
e = 0;
u = 0;
end;
vcov = [ 0.000081 0.0000081; 0.0000081 0.000081];
order = 1;
dynare++/tests/test2.mod deleted 100644 → 0
View file @ 219d2bb3
var y, c, k, a, h, b;
varexo e,u;
parameters beta, rho, alpha, delta, theta, psi, tau, phi;
alpha = 0.36;
rho = 0.95;
tau = 0.025;
beta = 0.99;
delta = 0.025;
psi = 0;
theta = 2.95;
phi = 0.1;
model;
c*theta*h^(1+psi)=(1-alpha)*y;
k = beta*(((exp(b)*c)/(exp(b(+1))*c(+1)))*(exp(b(+1))*alpha*y(+1)+(1-delta)*k));
y = exp(a)*(k(-1)^alpha)*(h^(1-alpha));
k = exp(b)*(y-c)+(1-delta)*k(-1);
a = rho*a(-1)+tau*b(-1) - rho*a(-2) - tau*b(-3) + e;
b = tau*a(-1)+rho*b(-1) - rho*b(-2) - tau*a(-3) + u;
end;
initval;
y = 1.08;
c = 0.8;
h = 0.29;
k = 11.08;
a = 0;
b = 0;
e = 0;
u = 0;
end;
vcov = [ 0.01 0.005; 0.005 0.01];
order = 1;
dynare++/tests/test2a.mod deleted 100644 → 0
View file @ 219d2bb3
var y, c, k, a, h, b;
varexo e,u;
parameters beta, rho, alpha, delta, theta, psi, tau, phi;
alpha = 0.36;
rho = 0.95;
tau = 0.025;
beta = 0.99;
delta = 0.025;
psi = 0;
theta = 2.95;
phi = 0.1;
model;
c*theta*h^(1+psi)=(1-alpha)*y;
k = beta*(((exp(b)*c)/(exp(b(+1))*c(+1)))*(exp(b(+1))*alpha*y(+1)+(1-delta)*k));
y = exp(a)*(k(-1)^alpha)*(h^(1-alpha));
k = exp(b)*(y-c)+(1-delta)*k(-1);
a = rho*a(-1)+tau*b(-1) - rho*a(-2) - tau*b(-3) + e;
b = tau*a(-1)+rho*b(-1) - rho*b(-2) - tau*a(-3) + u;
end;
initval;
y = 1.08;
c = 0.8;
h = 0.29;
k = 11.08;
a = 0;
b = 0;
e = 0;
u = 0;
end;
vcov = [ 0.01 0.005; 0.005 0.01];
order = 2;
dynare++/tests/test3.mod deleted 100644 → 0
View file @ 219d2bb3
var y,x;
varexo u,v;
parameters a, b, c, d, e, f, g, h, j;
a=0.8;
b=0.9;
c=0.9;
d=1;
e=-0.556875;
f=-0.172125;
g=-0.9;
h=-0.2754;
j=-1.8;
model;
x=a*x(-1)+u;
c*y(+1)^2+d*y^2+e*x^2+f*u^2-d*v^2+g+h*x(-1)*u+j*x(-1)*v=0;
end;
initval;
x=0;
y=0.7237469;
u=0;
v=0;
end;
vcov=[1 0; 0 1];
order = 2;
\ No newline at end of file
dynare++/tests/test4.mod deleted 100644 → 0
View file @ 219d2bb3
var y,x;
varexo u,v;
parameters a, b, c, d, e, f, m;
a=0.8;
b=0.9;
c=0.9;
d=1;
e=1;
m=50;
f = 1;
model;
x = a*x(-1)+u;
c*y(+1)^2+d*y^2+e*x^2-(c+d)*m^2-(c*b*b*a*a+d*b*b+e*a*a)*x(-1)^2-(c*b*b+e)*u^2-2*(c*m*b*a+d*m*b)*x(-1)-2*c*m*b*u-2*(c*b*b*a+e*a)*x(-1)*u-d*f^2*v^2-2*d*m*f*v-2*d*b*f*x(-1)*v=0;
end;
initval;
x=1;
y=21;
u=0;
v=0;
end;
vcov=[1 0; 0 1];
order = 2;
dynare++/tests/test5.mod deleted 100644 → 0
View file @ 219d2bb3
var y,x;
varexo u,v;
parameters a, b, c, d, e, m, n;
a=-0.8;
b=0.9;
c=0.9;
d=1;
e=1;
m=50;
n=0.2;
model;
x=b*x(-1)+u;
a*y(+1)+y-(a*b^2+1)*x(-1)^2-2*a*b*x(-1)*u-a*u^2-a-2*x(-1)*v-v^2;
end;
initval;
x=0;
y=0;
u=0;
v=0;
end;
vcov=[1 0; 0 1];
order = 3;
dynare++/tests/test6.mod deleted 100644 → 0
View file @ 219d2bb3
var y,x;
varexo u,v;
parameters a, b, c, d, e, m, n;
a=-0.8;
b=0.9;
c=0.9;
d=1;
e=1;
m=50;
n=0.2;
model;
x=b*x(-1)+u;
a*y(+1)+y-(a*b^3+1)*x(-1)^3-3*a*b*x(-1)*u^2-3*a*b^2*x(-1)^2*u-a*u^3-a-v^2;
end;
initval;
x=0;
y=0;
u=0;
v=0;
end;
vcov=[1 0; 0 1];
order = 3;
dynare++/tests/test7.mod deleted 100644 → 0
View file @ 219d2bb3
var y,x;
varexo u,v;
parameters a, b, c, d, e, m, n;
a=-0.8;
b=0.9;
c=0.9;
d=1;
e=1;
m=50;
n=0.2;
model;
x=b*x(-1)+u;
a*y(+1)+y-(a*b^4+1)*x(-1)^4-4*a*b*x(-1)*u^3-4*a*b^3*x(-1)^3*u-6*a*(b*x(-1)*u)^2-a*u^4-v;
end;
initval;
x=0;
y=0;
u=0;
v=0;
end;
vcov=[1 0; 0 1];
order = 4;
dynare++/tl/Makefile.am deleted 100644 → 0
View file @ 219d2bb3
SUBDIRS = cc testing
EXTRA_DIST = tl.tex
if HAVE_PDFLATEX
pdf-local: tl.pdf
tl.pdf: tl.tex
$(PDFLATEX) tl
$(PDFLATEX) tl
endif
CLEANFILES = *.pdf *.log *.aux
dynare++/tl/cc/Makefile.am deleted 100644 → 0
View file @ 219d2bb3
noinst_LIBRARIES = libtl.a
libtl_a_SOURCES = \
equivalence.cc \
equivalence.hh \
fine_container.cc \
fine_container.hh \
fs_tensor.cc \
fs_tensor.hh \
gs_tensor.cc \
gs_tensor.hh \
int_sequence.cc \
int_sequence.hh \
kron_prod.cc \
kron_prod.hh \
normal_moments.cc \
normal_moments.hh \
permutation.cc \
permutation.hh \
ps_tensor.cc \
ps_tensor.hh \
pyramid_prod.cc \
pyramid_prod.hh \
pyramid_prod2.cc \
pyramid_prod2.hh \
rfs_tensor.cc \
rfs_tensor.hh \
sparse_tensor.cc \
sparse_tensor.hh \
stack_container.cc \
stack_container.hh \
sthread.cc \
sthread.hh \
symmetry.cc \
symmetry.hh \
t_container.cc \
t_container.hh \
t_polynomial.cc \
t_polynomial.hh \
tensor.cc \
tensor.hh \
tl_exception.hh \
tl_static.cc \
tl_static.hh \
twod_matrix.cc \
twod_matrix.hh
libtl_a_CPPFLAGS = -I../../sylv/cc $(CPPFLAGS_MATIO)
libtl_a_CXXFLAGS = $(AM_CXXFLAGS) $(THREAD_CXXFLAGS)
dynare++/tl/cc/fs_tensor.hh deleted 100644 → 0
View file @ 219d2bb3
// Copyright 2004, Ondra Kamenik
// Full symmetry tensor.
/* Here we define folded and unfolded tensors for full symmetry. All
tensors from here are identifying the multidimensional index with
columns. */
#ifndef FS_TENSOR_H
#define FS_TENSOR_H
#include "tensor.hh"
#include "symmetry.hh"
class FGSTensor;
class UGSTensor;
class FRSingleTensor;
class FSSparseTensor;
/* Folded tensor with full symmetry maintains only information about
number of symmetrical variables |nv|. Further, we implement what is
left from the super class |FTensor|.
We implement |getOffset| which should be used with care since
its complexity.
We implement a method adding a given general symmetry tensor to the
full symmetry tensor supposing the variables of the general symmetry
tensor are stacked giving only one variable of the full symmetry
tensor. For instance, if $x=[y^T, u^T]^T$, then we can add tensor
$\left[g_{y^2u}\right]$ to tensor $g_{x^3}$. This is done in method
|addSubTensor|. Consult |@<|FGSTensor| class declaration@>| to know
what is general symmetry tensor. */
class UFSTensor;
class FFSTensor : public FTensor
{
int nv;
public:
/* Here are the constructors. The second constructor constructs a
tensor by one-dimensional contraction from the higher dimensional
tensor |t|. This is, it constructs a tensor
$$\left[g_{y^n}\right]_{\alpha_1\ldots\alpha_n}=
\left[t_{y^{n+1}}\right]_{\alpha_1\ldots\alpha_n\beta}[x]^\beta$$
See implementation |@<|FFSTensor| contraction constructor@>| for details.
The next constructor converts from sparse tensor (which is fully
symmetric and folded by nature).
The fourth constructs object from unfolded fully symmetric.
The fifth constructs a subtensor of selected rows. */
FFSTensor(int r, int nvar, int d)
: FTensor(along_col, IntSequence(d, nvar),
r, calcMaxOffset(nvar, d), d), nv(nvar)
{
}
FFSTensor(const FFSTensor &t, const ConstVector &x);
FFSTensor(const FSSparseTensor &t);
FFSTensor(const FFSTensor &ft)
= default;
FFSTensor(const UFSTensor &ut);
FFSTensor(int first_row, int num, FFSTensor &t)
: FTensor(first_row, num, t), nv(t.nv)
{
}
void increment(IntSequence &v) const override;
void decrement(IntSequence &v) const override;
UTensor&unfold() const override;
Symmetry
getSym() const
{
return Symmetry(dimen());
}
int getOffset(const IntSequence &v) const override;
void addSubTensor(const FGSTensor &t);
int
nvar() const
{
return nv;
}
static int calcMaxOffset(int nvar, int d);
};
/* Unfolded fully symmetric tensor is almost the same in structure as
|FFSTensor|, but the method |unfoldData|. It takes columns which also
exist in folded version and copies them to all their symmetrical
locations. This is useful when constructing unfolded tensor from
folded one. */
class UFSTensor : public UTensor
{
int nv;
public:
UFSTensor(int r, int nvar, int d)
: UTensor(along_col, IntSequence(d, nvar),
r, calcMaxOffset(nvar, d), d), nv(nvar)
{
}
UFSTensor(const UFSTensor &t, const ConstVector &x);
UFSTensor(const UFSTensor &ut)
= default;
UFSTensor(const FFSTensor &ft);
UFSTensor(int first_row, int num, UFSTensor &t)
: UTensor(first_row, num, t), nv(t.nv)
{
}
void increment(IntSequence &v) const override;
void decrement(IntSequence &v) const override;
FTensor&fold() const override;
Symmetry
getSym() const
{
return Symmetry(dimen());
}
int getOffset(const IntSequence &v) const override;
void addSubTensor(const UGSTensor &t);
int
nvar() const
{
return nv;
}
static int
calcMaxOffset(int nvar, int d)
{
return power(nvar, d);
}
private:
void unfoldData();
};
#endif
dynare++/tl/cc/int_sequence.cc deleted 100644 → 0
View file @ 219d2bb3
// Copyright 2004, Ondra Kamenik
#include "int_sequence.hh"
#include "symmetry.hh"
#include "tl_exception.hh"
#include <iostream>
#include <limits>
#include <numeric>
/* This unfolds a given integer sequence with respect to the given
symmetry. If for example the symmetry is $(2,3)$, and the sequence is
$(a,b)$, then the result is $(a,a,b,b,b)$. */
IntSequence::IntSequence(const Symmetry &sy, const IntSequence &se)
: data{new int[sy.dimen()], [](int *arr) { delete[] arr; }}, length{sy.dimen()}
{
int k = 0;
for (int i = 0; i < sy.num(); i++)
for (int j = 0; j < sy[i]; j++, k++)
operator[](k) = se[i];
}
/* This constructs an implied symmetry (implemented as |IntSequence|
from a more general symmetry and equivalence class (implemented as
|vector<int>|). For example, let the general symmetry be $y^3u^2$ and
the equivalence class is $\{0,4\}$ picking up first and fifth
variable, we calculate symmetry (at this point only |IntSequence|)
corresponding to the picked variables. These are $yu$. Thus the
constructed sequence must be $(1,1)$, meaning that we picked one $y$
and one $u$. */
IntSequence::IntSequence(const Symmetry &sy, const std::vector<int> &se)
: data{new int[sy.num()], [](int *arr) { delete[] arr; }}, length{sy.num()}
{
TL_RAISE_IF(sy.dimen() <= se[se.size()-1],
"Sequence is not reachable by symmetry in IntSequence()");
for (int i = 0; i < length; i++)
operator[](i) = 0;
for (int i : se)
operator[](sy.findClass(i))++;
}
/* This constructs an ordered integer sequence from the given ordered
sequence inserting the given number to the sequence. */
IntSequence::IntSequence(int i, const IntSequence &s)
: data{new int[s.size()+1], [](int *arr) { delete[] arr; }}, length{s.size()+1}
{
int j = 0;
while (j < s.size() && s[j] < i)
j++;
for (int jj = 0; jj < j; jj++)
operator[](jj) = s[jj];
operator[](j) = i;
for (int jj = j; jj < s.size(); jj++)
operator[](jj+1) = s[jj];
}
IntSequence::IntSequence(int i, const IntSequence &s, int pos)
: data{new int[s.size()+1], [](int *arr) { delete[] arr; }}, length{s.size()+1}
{
TL_RAISE_IF(pos < 0 || pos > s.size(),
"Wrong position for insertion IntSequence constructor");
for (int jj = 0; jj < pos; jj++)
operator[](jj) = s[jj];
operator[](pos) = i;
for (int jj = pos; jj < s.size(); jj++)
operator[](jj+1) = s[jj];
}
const IntSequence &
IntSequence::operator=(const IntSequence &s)
{
TL_RAISE_IF(length != s.length, "Wrong length for in-place IntSequence::operator=");
std::copy_n(s.data.get()+s.offset, length, data.get()+offset);
return *this;
}
const IntSequence &
IntSequence::operator=(IntSequence &&s)
{
TL_RAISE_IF(length != s.length, "Wrong length for in-place IntSequence::operator=");
std::copy_n(s.data.get()+s.offset, length, data.get()+offset);
return *this;
}
bool
IntSequence::operator==(const IntSequence &s) const
{
return std::equal(data.get()+offset, data.get()+offset+length,
s.data.get()+s.offset, s.data.get()+s.offset+s.length);
}
bool
IntSequence::operator<(const IntSequence &s) const
{
return std::lexicographical_compare(data.get()+offset, data.get()+offset+length,
s.data.get()+s.offset, s.data.get()+s.offset+s.length);
}
bool
IntSequence::lessEq(const IntSequence &s) const
{
TL_RAISE_IF(size() != s.size(),
"Sequence with different lengths in IntSequence::lessEq");
int i = 0;
while (i < size() && operator[](i) <= s[i])
i++;
return (i == size());
}
bool
IntSequence::less(const IntSequence &s) const
{
TL_RAISE_IF(size() != s.size(),
"Sequence with different lengths in IntSequence::less");
int i = 0;
while (i < size() && operator[](i) < s[i])
i++;
return (i == size());
}
void
IntSequence::sort()
{
std::sort(data.get()+offset, data.get()+offset+length);
}
/* Here we monotonize the sequence. If an item is less then its
predecessor, it is equalized. */
void
IntSequence::monotone()
{
for (int i = 1; i < length; i++)
if (operator[](i-1) > operator[](i))
operator[](i) = operator[](i-1);
}
/* This partially monotones the sequence. The partitioning is done by a
symmetry. So the subsequence given by the symmetry classes are
monotonized. For example, if the symmetry is $y^2u^3$, and the
|IntSequence| is $(5,3,1,6,4)$, the result is $(5,5,1,6,6)$. */
void
IntSequence::pmonotone(const Symmetry &s)
{
int cum = 0;
for (int i = 0; i < s.num(); i++)
{
for (int j = cum + 1; j < cum + s[i]; j++)
if (operator[](j-1) > operator[](j))
operator[](j) = operator[](j-1);
cum += s[i];
}
}
/* This returns sum of all elements. Useful for symmetries. */
int
IntSequence::sum() const
{
return std::accumulate(data.get()+offset, data.get()+offset+length, 0);
}
/* This returns product of subsequent items. Useful for Kronecker product
dimensions. */
int
IntSequence::mult(int i1, int i2) const
{
return std::accumulate(data.get()+offset+i1, data.get()+offset+i2,
1, std::multiplies<int>());
}
/* Return a number of the same items in the beginning of the sequence. */
int
IntSequence::getPrefixLength() const
{
int i = 0;
while (i+1 < size() && operator[](i+1) == operator[](0))
i++;
return i+1;
}
/* This returns a number of distinct items in the sequence. It supposes
that the sequence is ordered. For the empty sequence it returns zero. */
int
IntSequence::getNumDistinct() const
{
int res = 0;
if (length > 0)
res++;
for (int i = 1; i < length; i++)
if (operator[](i) != operator[](i-1))
res++;
return res;
}
/* This returns a maximum of the sequence. If the sequence is empty, it
returns the least possible |int| value. */
int
IntSequence::getMax() const
{
if (length == 0)
return std::numeric_limits<int>::min();
return *std::max_element(data.get()+offset, data.get()+offset+length);
}
void
IntSequence::add(int i)
{
for (int j = 0; j < size(); j++)
operator[](j) += i;
}
void
IntSequence::add(int f, const IntSequence &s)
{
TL_RAISE_IF(size() != s.size(),
"Wrong sequence length in IntSequence::add");
for (int j = 0; j < size(); j++)
operator[](j) += f*s[j];
}
bool
IntSequence::isPositive() const
{
return std::all_of(data.get()+offset, data.get()+offset+length,
[](int x) { return x >= 0; });
}
bool
IntSequence::isConstant() const
{
if (length < 2)
return true;
return std::all_of(data.get()+offset+1, data.get()+offset+length,
[this](int x) { return x == operator[](0); });
}
bool
IntSequence::isSorted() const
{
return std::is_sorted(data.get()+offset, data.get()+offset+length);
}
/* Debug print. */
void
IntSequence::print() const
{
std::cout << '[';
for (int i = 0; i < size(); i++)
std::cout << operator[](i) << ' ';
std::cout << ']' << std::endl;
}
dynare++/tl/cc/int_sequence.hh deleted 100644 → 0
View file @ 219d2bb3
// Copyright 2004, Ondra Kamenik
// Integer sequence.
/* Here we define an auxiliary abstraction for a sequence of integers. The
basic functionality is to hold an ordered sequence of integers with
constant length. We prefer using this simple class before STL
|vector<int>| since it is more efficient for our purposes.
The class is used in index of a tensor, in symmetry definition, in
Kronecker product dimensions, or as a class of an equivalence. The
latter case is not ordered, but we always order equivalence classes in
order to ensure unique representativeness. For almost all cases we
need the integer sequence to be ordered (sort), or monotonize (indices
of folded tensors), or partially monotonize (indices of folded tensors
not fully symmetric), or calculate a product of all members or only of
a part (used in Kronecker product dimensions). When we calculate
offsets in folded tensors, we need to obtain a number of the same
items in the front (|getPrefixLength|), and also to add some integer
number to all items.
Also, we need to construct a subsequence of a sequence. */
#ifndef INT_SEQUENCE_H
#define INT_SEQUENCE_H
#include <vector>
#include <memory>
#include <algorithm>
/* The implementation of |IntSequence| is straightforward. It has a pointer
|data|, an |offset| integer indicating the beginning of the data relatively
to the pointer and a |length| of the sequence. */
class Symmetry;
class IntSequence
{
std::shared_ptr<int> data;
int length;
int offset{0};
public:
/* We have a constructor allocating a given length of data, constructor
allocating and then initializing all members to a given number, a copy
constructor, a conversion from |vector<int>|, a subsequence
constructor, a constructor used for calculating implied symmetry from
a more general symmetry and one equivalence class (see |Symmetry|
class). Finally we have a constructor which unfolds a sequence with
respect to a given symmetry and constructor which inserts a given
number to the ordered sequence or given number to a given position. */
explicit IntSequence(int l)
: data{new int[l], [](int *arr) { delete[] arr; }}, length{l}
{
}
IntSequence(int l, int n)
: data{new int[l], [](int *arr) { delete[] arr; }}, length{l}
{
std::fill_n(data.get(), length, n);
}
IntSequence(const IntSequence &s)
: data{new int[s.length], [](int *arr) { delete[] arr; }}, length{s.length}
{
std::copy_n(s.data.get()+s.offset, length, data.get());
}
IntSequence(IntSequence &&s) = default;
IntSequence(IntSequence &s, int i1, int i2)
: data{s.data}, length{i2-i1}, offset{s.offset+i1}
{
}
IntSequence(const IntSequence &s, int i1, int i2)
: data{new int[i2-i1], [](int *arr) { delete[] arr; }}, length{i2-i1}
{
std::copy_n(s.data.get()+s.offset+i1, length, data.get());
}
IntSequence(const Symmetry &sy, const std::vector<int> &se);
IntSequence(const Symmetry &sy, const IntSequence &se);
IntSequence(int i, const IntSequence &s);
IntSequence(int i, const IntSequence &s, int pos);
const IntSequence &operator=(const IntSequence &s);
const IntSequence &operator=(IntSequence &&s);
virtual ~IntSequence() = default;
bool operator==(const IntSequence &s) const;
bool
operator!=(const IntSequence &s) const
{
return !operator==(s);
}
int &
operator[](int i)
{
return data.get()[offset+i];
}
int
operator[](int i) const
{
return data.get()[offset+i];
}
int
size() const
{
return length;
}
/* We provide two orderings. The first |operator<| is the linear
lexicographic ordering, the second |less| is the non-linear Cartesian
ordering. */
bool operator<(const IntSequence &s) const;
bool
operator<=(const IntSequence &s) const
{
return (operator==(s) || operator<(s));
}
bool lessEq(const IntSequence &s) const;
bool less(const IntSequence &s) const;
void sort();
void monotone();
void pmonotone(const Symmetry &s);
int sum() const;
int mult(int i1, int i2) const;
int
mult() const
{
return mult(0, length);
}
void add(int i);
void add(int f, const IntSequence &s);
int getPrefixLength() const;
int getNumDistinct() const;
int getMax() const;
bool isPositive() const;
bool isConstant() const;
bool isSorted() const;
void print() const;
};
#endif
dynare++/tl/cc/kron_prod.cc deleted 100644 → 0
View file @ 219d2bb3
// Copyright 2004, Ondra Kamenik
#include "kron_prod.hh"
#include "tl_exception.hh"
#include <cstdio>
/* Here we construct Kronecker product dimensions from Kronecker
product dimensions by picking a given matrix and all other set to
identity. The constructor takes dimensions of $A_1\otimes
A_2\otimes\ldots\otimes A_n$, and makes dimensions of $I\otimes
A_i\otimes I$, or $I\otimes A_n$, or $A_1\otimes I$ for a given
$i$. The identity matrices must fit into the described order. See
header file.
We first decide what is a length of the resulting dimensions. Possible
length is three for $I\otimes A\otimes I$, and two for $I\otimes A$,
or $A\otimes I$.
Then we fork according to |i|. */
KronProdDimens::KronProdDimens(const KronProdDimens &kd, int i)
: rows((i == 0 || i == kd.dimen()-1) ? (2) : (3)),
cols((i == 0 || i == kd.dimen()-1) ? (2) : (3))
{
TL_RAISE_IF(i < 0 || i >= kd.dimen(),
"Wrong index for pickup in KronProdDimens constructor");
int kdim = kd.dimen();
if (i == 0)
{
// set AI dimensions
/* The first rows and cols are taken from |kd|. The dimensions of
identity matrix is a number of rows in $A_2\otimes\ldots\otimes A_n$
since the matrix $A_1\otimes I$ is the first. */
rows[0] = kd.rows[0];
rows[1] = kd.rows.mult(1, kdim);
cols[0] = kd.cols[0];
cols[1] = rows[1];
}
else if (i == kdim-1)
{
// set IA dimensions
/* The second dimension is taken from |kd|. The dimensions of identity
matrix is a number of columns of $A_1\otimes\ldots A_{n-1}$, since the
matrix $I\otimes A_n$ is the last. */
rows[0] = kd.cols.mult(0, kdim-1);
rows[1] = kd.rows[kdim-1];
cols[0] = rows[0];
cols[1] = kd.cols[kdim-1];
}
else
{
// set IAI dimensions
/* The dimensions of the middle matrix are taken from |kd|. The
dimensions of the first identity matrix are a number of columns of
$A_1\otimes\ldots\otimes A_{i-1}$, and the dimensions of the last
identity matrix are a number of rows of $A_{i+1}\otimes\ldots\otimes
A_n$. */
rows[0] = kd.cols.mult(0, i);
cols[0] = rows[0];
rows[1] = kd.rows[i];
cols[1] = kd.cols[i];
cols[2] = kd.rows.mult(i+1, kdim);
rows[2] = cols[2];
}
}
/* This raises an exception if dimensions are bad for multiplication
|out = in*this|. */
void
KronProd::checkDimForMult(const ConstTwoDMatrix &in, const TwoDMatrix &out) const
{
int my_rows;
int my_cols;
kpd.getRC(my_rows, my_cols);
TL_RAISE_IF(in.nrows() != out.nrows() || in.ncols() != my_rows,
"Wrong dimensions for KronProd in KronProd::checkDimForMult");
}
/* Here we Kronecker multiply two given vectors |v1| and |v2| and
store the result in preallocated |res|. */
void
KronProd::kronMult(const ConstVector &v1, const ConstVector &v2,
Vector &res)
{
TL_RAISE_IF(res.length() != v1.length()*v2.length(),
"Wrong vector lengths in KronProd::kronMult");
res.zeros();
for (int i = 0; i < v1.length(); i++)
{
Vector sub(res, i *v2.length(), v2.length());
sub.add(v1[i], v2);
}
}
void
KronProdAll::setMat(int i, const TwoDMatrix &m)
{
matlist[i] = &m;
kpd.setRC(i, m.nrows(), m.ncols());
}
void
KronProdAll::setUnit(int i, int n)
{
matlist[i] = nullptr;
kpd.setRC(i, n, n);
}
bool
KronProdAll::isUnit() const
{
int i = 0;
while (i < dimen() && matlist[i] == nullptr)
i++;
return i == dimen();
}
/* Here we multiply $B\cdot(I\otimes A)$. If $m$ is a dimension of the
identity matrix, then the product is equal to
$B\cdot\hbox{diag}_m(A)$. If $B$ is partitioned accordingly, then the
result is $[B_1A, B_2A,\ldots B_mA]$.
Here, |outi| are partitions of |out|, |ini| are const partitions of
|in|, and |id_cols| is $m$. We employ level-2 BLAS. */
void
KronProdIA::mult(const ConstTwoDMatrix &in, TwoDMatrix &out) const
{
checkDimForMult(in, out);
int id_cols = kpd.cols[0];
ConstTwoDMatrix a(mat);
for (int i = 0; i < id_cols; i++)
{
TwoDMatrix outi(out, i *a.ncols(), a.ncols());
ConstTwoDMatrix ini(in, i *a.nrows(), a.nrows());
outi.mult(ini, a);
}
}
/* Here we construct |KronProdAI| from |KronProdIAI|. It is clear. */
KronProdAI::KronProdAI(const KronProdIAI &kpiai)
: KronProd(KronProdDimens(2)), mat(kpiai.mat)
{
kpd.rows[0] = mat.nrows();
kpd.cols[0] = mat.ncols();
kpd.rows[1] = kpiai.kpd.rows[2];
kpd.cols[1] = kpiai.kpd.cols[2];
}
/* Here we multiply $B\cdot(A\otimes I)$. Let the dimension of the
matrix $A$ be $m\times n$, the dimension of $I$ be $p$, and a number
of rows of $B$ be $q$. We use the fact that $B\cdot(A\otimes
I)=\hbox{reshape}(\hbox{reshape}(B, q, mp)\cdot A, q, np)$. This works
only for matrix $B$, whose storage has leading dimension equal to
number of rows.
For cases where the leading dimension is not equal to the number of
rows, we partition the matrix $A\otimes I$ to $m\times n$ square
partitions $a_{ij}I$. Therefore, we partition $B$ to $m$ partitions
$[B_1, B_2,\ldots,B_m]$. Each partition of $B$ has the same number of
columns as the identity matrix. If $R$ denotes the resulting matrix,
then it can be partitioned to $n$ partitions
$[R_1,R_2,\ldots,R_n]$. Each partition of $R$ has the same number of
columns as the identity matrix. Then we have $R_i=\sum a_{ji}B_j$.
In code, |outi| is $R_i$, |ini| is $B_j$, and |id_cols| is a dimension
of the identity matrix */
void
KronProdAI::mult(const ConstTwoDMatrix &in, TwoDMatrix &out) const
{
checkDimForMult(in, out);
int id_cols = kpd.cols[1];
ConstTwoDMatrix a(mat);
if (in.getLD() == in.nrows())
{
ConstTwoDMatrix in_resh(in.nrows()*id_cols, a.nrows(), in.getData());
TwoDMatrix out_resh(in.nrows()*id_cols, a.ncols(), out.getData());
out_resh.mult(in_resh, a);
}
else
{
out.zeros();
for (int i = 0; i < a.ncols(); i++)
{
TwoDMatrix outi(out, i *id_cols, id_cols);
for (int j = 0; j < a.nrows(); j++)
{
ConstTwoDMatrix ini(in, j *id_cols, id_cols);
outi.add(a.get(j, i), ini);
}
}
}
}
/* Here we multiply $B\cdot(I\otimes A\otimes I)$. If $n$ is a
dimension of the first identity matrix, then we multiply
$B\cdot\hbox{diag}_n(A\otimes I)$. So we partition $B$ and result $R$
accordingly, and multiply $B_i\cdot(A\otimes I)$, which is in fact
|KronProdAI::mult|. Note that number of columns of partitions of $B$
are number of rows of $A\otimes I$, and number of columns of $R$ are
number of columns of $A\otimes I$.
In code, |id_cols| is $n$, |akronid| is a Kronecker product object of
$A\otimes I$, and |in_bl_width|, and |out_bl_width| are rows and cols of
$A\otimes I$. */
void
KronProdIAI::mult(const ConstTwoDMatrix &in, TwoDMatrix &out) const
{
checkDimForMult(in, out);
int id_cols = kpd.cols[0];
KronProdAI akronid(*this);
int in_bl_width;
int out_bl_width;
akronid.kpd.getRC(in_bl_width, out_bl_width);
for (int i = 0; i < id_cols; i++)
{
TwoDMatrix outi(out, i *out_bl_width, out_bl_width);
ConstTwoDMatrix ini(in, i *in_bl_width, in_bl_width);
akronid.mult(ini, outi);
}
}
/* Here we multiply $B\cdot(A_1\otimes\ldots\otimes A_n)$. First we
multiply $B\cdot(A_1\otimes I)$, then this is multiplied by all
$I\otimes A_i\otimes I$, and finally by $I\otimes A_n$.
If the dimension of the Kronecker product is only 1, then we multiply
two matrices in straight way and return.
The intermediate results are stored on heap pointed by |last|. A new
result is allocated, and then the former storage is deallocated.
We have to be careful in cases when last or first matrix is unit and
no calculations are performed in corresponding codes. The codes should
handle |last| safely also if no calcs are done. */
void
KronProdAll::mult(const ConstTwoDMatrix &in, TwoDMatrix &out) const
{
// quick copy if product is unit
if (isUnit())
{
out.zeros();
out.add(1.0, in);
return;
}
// quick zero if one of the matrices is zero
/* If one of the matrices is exactly zero or the |in| matrix is zero,
set out to zero and return */
bool is_zero = false;
for (int i = 0; i < dimen() && !is_zero; i++)
is_zero = matlist[i] && matlist[i]->isZero();
if (is_zero || in.isZero())
{
out.zeros();
return;
}
// quick multiplication if dimension is 1
if (dimen() == 1)
{
if (matlist[0]) // always true
out.mult(in, ConstTwoDMatrix(*(matlist[0])));
return;
}
int c;
TwoDMatrix *last = nullptr;
// perform first multiplication AI
/* Here we have to construct $A_1\otimes I$, allocate intermediate
result |last|, and perform the multiplication. */
if (matlist[0])
{
KronProdAI akronid(*this);
c = akronid.kpd.ncols();
last = new TwoDMatrix(in.nrows(), c);
akronid.mult(in, *last);
}
else
{
last = new TwoDMatrix(in.nrows(), in.ncols(), Vector{in.getData()});
}
// perform intermediate multiplications IAI
/* Here we go through all $I\otimes A_i\otimes I$, construct the
product, allocate new storage for result |newlast|, perform the
multiplication, deallocate old |last|, and set |last| to |newlast|. */
for (int i = 1; i < dimen()-1; i++)
{
if (matlist[i])
{
KronProdIAI interkron(*this, i);
c = interkron.kpd.ncols();
auto *newlast = new TwoDMatrix(in.nrows(), c);
interkron.mult(*last, *newlast);
delete last;
last = newlast;
}
}
// perform last multiplication IA
/* Here just construct $I\otimes A_n$ and perform multiplication and
deallocate |last|. */
if (matlist[dimen()-1])
{
KronProdIA idkrona(*this);
idkrona.mult(*last, out);
}
else
{
out = *last;
}
delete last;
}
/* This calculates a Kornecker product of rows of matrices, the row
indices are given by the integer sequence. The result is allocated and
returned. The caller is repsonsible for its deallocation. */
Vector *
KronProdAll::multRows(const IntSequence &irows) const
{
TL_RAISE_IF(irows.size() != dimen(),
"Wrong length of row indices in KronProdAll::multRows");
Vector *last = nullptr;
ConstVector *row;
std::vector<Vector *> to_delete;
for (int i = 0; i < dimen(); i++)
{
int j = dimen()-1-i;
// set |row| to the row of |j|-th matrix
/* If the |j|-th matrix is real matrix, then the row is constructed
from the matrix. It the matrix is unit, we construct a new vector,
fill it with zeros, than set the unit to appropriate place, and make
the |row| as ConstVector of this vector, which sheduled for
deallocation. */
if (matlist[j])
row = new ConstVector(matlist[j]->getRow(irows[j]));
else
{
auto *aux = new Vector(ncols(j));
aux->zeros();
(*aux)[irows[j]] = 1.0;
to_delete.push_back(aux);
row = new ConstVector(*aux);
}
// set |last| to product of |row| and |last|
/* If the |last| is exists, we allocate new storage, Kronecker
multiply, deallocate the old storage. If the |last| does not exist,
then we only make |last| equal to |row|. */
if (last)
{
Vector *newlast;
newlast = new Vector(last->length()*row->length());
kronMult(*row, ConstVector(*last), *newlast);
delete last;
last = newlast;
}
else
{
last = new Vector(*row);
}
delete row;
}
for (auto & i : to_delete)
delete i;
return last;
}
/* This permutes the matrices so that the new ordering would minimize
memory consumption. As shown in |@<|KronProdAllOptim| class declaration@>|,
we want ${m_k\over n_k}\leq{m_{k-1}\over n_{k-1}}\ldots\leq{m_1\over n_1}$,
where $(m_i,n_i)$ is the dimension of $A_i$. So we implement the bubble
sort. */
void
KronProdAllOptim::optimizeOrder()
{
for (int i = 0; i < dimen(); i++)
{
int swaps = 0;
for (int j = 0; j < dimen()-1; j++)
{
if (((double) kpd.rows[j])/kpd.cols[j] < ((double) kpd.rows[j+1])/kpd.cols[j+1])
{
// swap dimensions and matrices at |j| and |j+1|
int s = kpd.rows[j+1];
kpd.rows[j+1] = kpd.rows[j];
kpd.rows[j] = s;
s = kpd.cols[j+1];
kpd.cols[j+1] = kpd.cols[j];
kpd.cols[j] = s;
const TwoDMatrix *m = matlist[j+1];
matlist[j+1] = matlist[j];
matlist[j] = m;
// project the swap to the permutation |oper|
s = oper.getMap()[j+1];
oper.getMap()[j+1] = oper.getMap()[j];
oper.getMap()[j] = s;
swaps++;
}
}
if (swaps == 0)
{
return;
}
}
}
dynare++/tl/cc/kron_prod.hh deleted 100644 → 0
View file @ 219d2bb3
// Copyright 2004, Ondra Kamenik
// Kronecker product.
/* Here we define an abstraction for a Kronecker product of a sequence of
matrices. This is $A_1\otimes\ldots\otimes A_n$. Obviously we do not
store the product in memory. First we need to represent a dimension
of the Kronecker product. Then we represent the Kronecker product,
simply it is the Kronecker product dimension with a vector of
references to the matrices $A_1,\ldots, A_n$.
The main task of this class is to calculate a matrix product
$B\cdot(A_1\otimes A_2\otimes\ldots\otimes A_n)$ which in
our application has much more moderate dimensions than $A_1\otimes
A_2\otimes\ldots\otimes A_n$. We calculate it as
$$B\cdot(A_1\otimes I)\cdot\ldots\cdot(I\otimes A_i\otimes
I)\cdot\ldots\cdot (I\otimes A_n)$$
where dimensions of identity matrices differ and are given by the
chosen order. One can naturally ask, whether there is some optimal
order minimizing maximum storage needed for intermediate
results. The optimal ordering is implemented by class |KronProdAllOptim|.
For this multiplication, we also need to represent products of type
$A\otimes I$, $I\otimes A\otimes I$, and $I\otimes A$. */
#ifndef KRON_PROD_H
#define KRON_PROD_H
#include "twod_matrix.hh"
#include "permutation.hh"
#include "int_sequence.hh"
class KronProdAll;
class KronProdAllOptim;
class KronProdIA;
class KronProdIAI;
class KronProdAI;
/* |KronProdDimens| maintains a dimension of the Kronecker product. So,
it maintains two sequences, one for rows, and one for columns. */
class KronProdDimens
{
friend class KronProdAll;
friend class KronProdAllOptim;
friend class KronProdIA;
friend class KronProdIAI;
friend class KronProdAI;
private:
IntSequence rows;
IntSequence cols;
public:
/* We define three constructors. First initializes to a given
dimension, and all rows and cols are set to zeros. Second is a copy
constructor. The third constructor takes dimensions of $A_1\otimes
A_2\otimes\ldots\otimes A_n$, and makes dimensions of $I\otimes
A_i\otimes I$, or $I\otimes A_n$, or $A_1\otimes I$ for a given
$i$. The dimensions of identity matrices are such that
$$A_1\otimes A_2\otimes\ldots\otimes A_n=
(A_1\otimes I)\cdot\ldots\cdot(I\otimes A_i\otimes I)
\cdot\ldots\cdot(I\otimes A_n)$$
Note that the matrices on the right do not commute only because sizes
of identity matrices which are then given by this ordering. */
KronProdDimens(int dim)
: rows(dim, 0), cols(dim, 0)
{
}
KronProdDimens(const KronProdDimens &kd)
= default;
KronProdDimens(const KronProdDimens &kd, int i);
KronProdDimens &
operator=(const KronProdDimens &kd)
= default;
bool
operator==(const KronProdDimens &kd) const
{
return rows == kd.rows && cols == kd.cols;
}
int
dimen() const
{
return rows.size();
}
void
setRC(int i, int r, int c)
{
rows[i] = r; cols[i] = c;
}
void
getRC(int i, int &r, int &c) const
{
r = rows[i]; c = cols[i];
}
void
getRC(int &r, int &c) const
{
r = rows.mult(); c = cols.mult();
}
int
nrows() const
{
return rows.mult();
}
int
ncols() const
{
return cols.mult();
}
int
nrows(int i) const
{
return rows[i];
}
int
ncols(int i) const
{
return cols[i];
}
};
/* Here we define an abstract class for all Kronecker product classes,
which are |KronProdAll| (the most general), |KronProdIA| (for
$I\otimes A$), |KronProdAI| (for $A\otimes I$), and |KronProdIAI| (for
$I\otimes A\otimes I$). The purpose of the super class is to only
define some common methods and common member |kpd| for dimensions and
declare pure virtual |mult| which is implemented by the subclasses.
The class also contains a static method |kronMult|, which calculates a
Kronecker product of two vectors and stores it in the provided
vector. It is useful at a few points of the library. */
class KronProd
{
protected:
KronProdDimens kpd;
public:
KronProd(int dim)
: kpd(dim)
{
}
KronProd(const KronProdDimens &kd)
: kpd(kd)
{
}
KronProd(const KronProd &kp)
= default;
virtual ~KronProd()
= default;
int
dimen() const
{
return kpd.dimen();
}
virtual void mult(const ConstTwoDMatrix &in, TwoDMatrix &out) const = 0;
void
mult(const TwoDMatrix &in, TwoDMatrix &out) const
{
mult(ConstTwoDMatrix(in), out);
}
void checkDimForMult(const ConstTwoDMatrix &in, const TwoDMatrix &out) const;
void
checkDimForMult(const TwoDMatrix &in, const TwoDMatrix &out) const
{
checkDimForMult(ConstTwoDMatrix(in), out);
}
static void kronMult(const ConstVector &v1, const ConstVector &v2,
Vector &res);
int
nrows() const
{
return kpd.nrows();
}
int
ncols() const
{
return kpd.ncols();
}
int
nrows(int i) const
{
return kpd.nrows(i);
}
int
ncols(int i) const
{
return kpd.ncols(i);
}
};
/* |KronProdAll| is a main class of this file. It represents the
Kronecker product $A_1\otimes A_2\otimes\ldots\otimes A_n$. Besides
dimensions, it stores pointers to matrices in |matlist| array. If a
pointer is null, then the matrix is considered to be unit. The array
is set by calls to |setMat| method (for real matrices) or |setUnit|
method (for unit matrices).
The object is constructed by a constructor, which allocates the
|matlist| and initializes dimensions to zeros. Then a caller must feed
the object with matrices by calling |setMat| and |setUnit| repeatedly
for different indices.
We implement the |mult| method of |KronProd|, and a new method
|multRows|, which creates a vector of kronecker product of all rows of
matrices in the object. The rows are given by the |IntSequence|. */
class KronProdAll : public KronProd
{
friend class KronProdIA;
friend class KronProdIAI;
friend class KronProdAI;
protected:
const TwoDMatrix **const matlist;
public:
KronProdAll(int dim)
: KronProd(dim), matlist(new const TwoDMatrix *[dim])
{
}
~KronProdAll() override
{
delete [] matlist;
}
void setMat(int i, const TwoDMatrix &m);
void setUnit(int i, int n);
const TwoDMatrix &
getMat(int i) const
{
return *(matlist[i]);
}
void mult(const ConstTwoDMatrix &in, TwoDMatrix &out) const override;
Vector *multRows(const IntSequence &irows) const;
private:
bool isUnit() const;
};
/* The class |KronProdAllOptim| minimizes memory consumption of the
product $B\cdot(A_1\otimes A_2\otimes\ldots\otimes A_k)$. The
optimization is done by reordering of the matrices $A_1,\ldots,A_k$,
in order to minimize a sum of all storages needed for intermediate
results. The optimal ordering is also nearly optimal with respect to
number of flops.
Let $(m_i,n_i)$ be dimensions of $A_i$. It is easy to observe, that
for $i$-th step we need storage of $r\cdot n_1\cdot\ldots\cdot
n_i\cdot m_{i+1}\cdot\ldots\cdot m_k$, where $r$ is a number of rows
of $B$. To minimize the sum through all $i$ over all permutations of
matrices, it is equivalent to minimize the sum
$\sum_{i=1}^k{m_{i+1}\cdot\ldots\cdot m_k\over n_{i+1}\cdot\ldots\cdot
n_k}$. The optimal ordering will yield ${m_k\over
n_k}\leq{m_{k-1}\over n_{k-1}}\ldots\leq{m_1\over n_1}$.
Now observe, that the number of flops for $i$-th step is $r\cdot
n_1\cdot\ldots\cdot n_i\cdot m_i\cdot\ldots\cdot m_k$. In order to
minimize a number of flops, it is equivalent to minimize
$\sum_{i=1}^km_i{m_{i+1}\cdot\ldots\cdot m_k\over
n_{i+1}\cdot\ldots\cdot n_k}$. Note that, normally, the $m_i$ does not
change as much as $n_{j+1},\ldots,n_k$, so the ordering minimizing the
memory will be nearly optimal with respect to number of flops.
The class |KronProdAllOptim| inherits from |KronProdAll|. A public
method |optimizeOrder| does the reordering. The permutation is stored
in |oper|. So, as long as |optimizeOrder| is not called, the class is
equivalent to |KronProdAll|. */
class KronProdAllOptim : public KronProdAll
{
protected:
Permutation oper;
public:
KronProdAllOptim(int dim)
: KronProdAll(dim), oper(dim)
{
}
void optimizeOrder();
const Permutation &
getPer() const
{
return oper;
}
};
/* This class represents $I\otimes A$. We have only one reference to
the matrix, which is set by constructor. */
class KronProdIA : public KronProd
{
friend class KronProdAll;
const TwoDMatrix &mat;
public:
KronProdIA(const KronProdAll &kpa)
: KronProd(KronProdDimens(kpa.kpd, kpa.dimen()-1)),
mat(kpa.getMat(kpa.dimen()-1))
{
}
void mult(const ConstTwoDMatrix &in, TwoDMatrix &out) const override;
};
/* This class represents $A\otimes I$. We have only one reference to
the matrix, which is set by constructor. */
class KronProdAI : public KronProd
{
friend class KronProdIAI;
friend class KronProdAll;
const TwoDMatrix &mat;
public:
KronProdAI(const KronProdAll &kpa)
: KronProd(KronProdDimens(kpa.kpd, 0)),
mat(kpa.getMat(0))
{
}
KronProdAI(const KronProdIAI &kpiai);
void mult(const ConstTwoDMatrix &in, TwoDMatrix &out) const override;
};
/* This class represents $I\otimes A\otimes I$. We have only one reference to
the matrix, which is set by constructor. */
class KronProdIAI : public KronProd
{
friend class KronProdAI;
friend class KronProdAll;
const TwoDMatrix &mat;
public:
KronProdIAI(const KronProdAll &kpa, int i)
: KronProd(KronProdDimens(kpa.kpd, i)),
mat(kpa.getMat(i))
{
}
void mult(const ConstTwoDMatrix &in, TwoDMatrix &out) const override;
};
#endif
dynare++/tl/cc/normal_moments.cc deleted 100644 → 0
View file @ 219d2bb3
// Copyright 2004, Ondra Kamenik
#include "normal_moments.hh"
#include "permutation.hh"
#include "kron_prod.hh"
#include "tl_static.hh"
UNormalMoments::UNormalMoments(int maxdim, const TwoDMatrix &v)
: TensorContainer<URSingleTensor>(1)
{
if (maxdim >= 2)
generateMoments(maxdim, v);
}
/* Here we fill up the container with the tensors for $d=2,4,6,\ldots$
up to the given dimension. Each tensor of moments is equal to
$F_n\left(\otimes^nv\right).$ This has a dimension equal to
$2n$. See the header file for proof and details.
Here we sequentially construct the Kronecker power
$\otimes^nv$, and apply $F_n$. */
void
UNormalMoments::generateMoments(int maxdim, const TwoDMatrix &v)
{
TL_RAISE_IF(v.nrows() != v.ncols(),
"Variance-covariance matrix is not square in UNormalMoments constructor");
int nv = v.nrows();
auto *mom2 = new URSingleTensor(nv, 2);
mom2->getData() = v.getData();
insert(mom2);
auto *kronv = new URSingleTensor(nv, 2);
kronv->getData() = v.getData();
for (int d = 4; d <= maxdim; d += 2)
{
auto *newkronv = new URSingleTensor(nv, d);
KronProd::kronMult(ConstVector(v.getData()),
ConstVector(kronv->getData()),
newkronv->getData());
delete kronv;
kronv = newkronv;
auto *mom = new URSingleTensor(nv, d);
// apply $F_n$ to |kronv|
/* Here we go through all equivalences, select only those having 2
elements in each class, then go through all elements in |kronv| and
add to permuted location of |mom|.
The permutation must be taken as inverse of the permutation implied by
the equivalence, since we need a permutation which after application
to identity of indices yileds indices in the equivalence classes. Note
how the |Equivalence::apply| method works. */
mom->zeros();
const EquivalenceSet eset = ebundle.get(d);
for (const auto & cit : eset)
{
if (selectEquiv(cit))
{
Permutation per(cit);
per.inverse();
for (Tensor::index it = kronv->begin(); it != kronv->end(); ++it)
{
IntSequence ind(kronv->dimen());
per.apply(it.getCoor(), ind);
Tensor::index it2(mom, ind);
mom->get(*it2, 0) += kronv->get(*it, 0);
}
}
}
insert(mom);
}
delete kronv;
}
/* We return |true| for an equivalence whose each class has 2 elements. */
bool
UNormalMoments::selectEquiv(const Equivalence &e)
{
if (2*e.numClasses() != e.getN())
return false;
for (const auto & si : e)
{
if (si.length() != 2)
return false;
}
return true;
}
/* Here we go through all the unfolded container, fold each tensor and
insert it. */
FNormalMoments::FNormalMoments(const UNormalMoments &moms)
: TensorContainer<FRSingleTensor>(1)
{
for (const auto & mom : moms)
{
auto *fm = new FRSingleTensor(*(mom.second));
insert(fm);
}
}
dynare++/tl/cc/normal_moments.hh deleted 100644 → 0
View file @ 219d2bb3
// Copyright 2004, Ondra Kamenik
// Moments of normal distribution.
/* Here we calculate the higher order moments of normally distributed
random vector $u$ with means equal to zero and given
variance--covariance matrix $V$, this is $u\sim N(0,V)$. The moment
generating function for such distribution is $f(t)=e^{{1\over 2}t^TVt}$. If
we derivate it wrt $t$ and unfold the higher dimensional tensors
row-wise, we obtain terms like
$$\eqalign{
{\partial\over\partial t}f(t)=&f(t)\cdot Vt\cr
{\partial^2\over\partial t^2}f(t)=&f(t)\cdot(Vt\otimes Vt+v)\cr
{\partial^3\over\partial t^3}f(t)=&f(t)\cdot
(Vt\otimes Vt\otimes Vt+P_?(v\otimes Vt)+P_?(Vt\otimes v)+v\otimes Vt)\cr
{\partial^4\over\partial t^4}f(t)=&f(t)\cdot
(Vt\otimes Vt\otimes Vt\otimes Vt+S_?(v\otimes Vt\otimes Vt)+
S_?(Vt\otimes v\otimes Vt)+S_?(Vt\otimes Vt\otimes v)+S_?(v\otimes v))}
$$
where $v$ is vectorized $V$ ($v=\hbox{vec}(V)$), and $P_?$ is a
suitable row permutation (corresponds to permutation of
multidimensional indices) which permutes the tensor data, so that the
index of a variable being derived would be the last. This ensures that
all (permuted) tensors can be summed yielding a tensor whose indices
have some order (in here we chose the order that more recent
derivating variables are to the right). Finally, $S_?$ is a suitable
sum of various $P_?$.
We are interested in $S_?$ multiplying the Kronecker powers
$\otimes^nv$. The $S_?$ is a (possibly) multi-set of permutations of
even order. Note that we know a number of permutations in $S_?$. The
above formulas for $F(t)$ derivatives are valid also for monomial
$u$, and from literature we know that $2n$-th moment is ${(2n!)\over
n!2^n}\sigma^2$. So there are ${(2n!)\over n!2^n}$ permutations in
$S_?$.
In order to find the $S_?$ we need to define a couple of
things. First we define a sort of equivalence between the permutations
applicable to even number of indices. We write $P_1\equiv P_2$
whenever $P_1^{-1}\circ P_2$ permutes only whole pairs, or items
within pairs, but not indices across the pairs. For instance the
permutations $(0,1,2,3)$ and $(3,2,0,1)$ are equivalent, but
$(0,2,1,3)$ is not equivalent with the two. Clearly, the $\equiv$ is
an equivalence.
This allows to define a relation $\sqsubseteq$ between the permutation
multi-sets $S$, which is basically the subset relation $\subseteq$ but
with respect to the equivalence $\equiv$, more formally:
$$S_1\sqsubseteq S_2\quad\hbox{iff}\quad P\in S_1
\Rightarrow\exists Q\in S_2:P\equiv Q$$
This induces an equivalence $S_1\equiv S_2$.
Now let $F_n$ denote a set of permutations on $2n$ indices which is
maximal with respect to $\sqsubseteq$, and minimal with respect to
$\equiv$. (In other words, it contains everything up to the
equivalence $\equiv$.) It is straightforward to calculate a number of
permutations in $F_n$. This is a total number of all permutations of
$2n$ divided by permutations of pairs divided by permutations within
the pairs. This is ${(2n!)\over n!2^n}$.
We prove that $S_?\equiv F_n$. Clearly $S_?\sqsubseteq F_n$, since
$F_n$ is maximal. In order to prove that $F_n\sqsubseteq S_?$, let us
assert that for any permutation $P$ and for any (semi)positive
definite matrix $V$ we have $PS_?\otimes^nv=S_?\otimes^nv$. Below we
show that there is a positive definite matrix $V$ of some dimension
that for any two permutation multi-sets $S_1$, $S_2$, we have
$$S_1\not\equiv S_2\Rightarrow S_1(\otimes^nv)\neq S_2(\otimes^nv)$$
So it follows that for any permutation $P$, we have $PS_?\equiv
S_?$. For a purpose of contradiction let $P\in F_n$ be a permutation
which is not equivalent to any permutation from $S_?$. Since $S_?$ is
non-empty, let us pick $P_0\in S_?$. Now assert that
$P_0^{-1}S_?\not\equiv P^{-1}S_?$ since the first contains an identity
and the second does not contain a permutation equivalent to
identity. Thus we have $(P\circ P_0^{-1})S_?\not\equiv S_?$ which
gives the contradiction and we have proved that $F_n\sqsubseteq
S_?$. Thus $F_n\equiv S_?$. Moreover, we know that $S_?$ and $F_n$
have the same number of permutations, hence the minimality of $S_?$
with respect to $\equiv$.
Now it suffices to prove that there exists a positive definite $V$
such that for any two permutation multi-sets $S_1$, and $S_2$ holds
$S_1\not\equiv S_2\Rightarrow S_1(\otimes^nv)\neq S_2(\otimes^nv)$. If
$V$ is $n\times n$ matrix, then $S_1\not\equiv S_2$ implies that there
is identically nonzero polynomial of elements from $V$ of order $n$
over integers. If $V=A^TA$ then there is identically non-zero
polynomial of elements from $A$ of order $2n$. This means, that we
have to find $n(n+1)/2$ tuple $x$ of real numbers such that all
identically non-zero polynomials $p$ of order $2n$ over integers yield
$p(x)\neq 0$.
The $x$ is constructed as follows: $x_i = \pi^{\log{r_i}}$, where $r_i$
is $i$-th prime. Let us consider monom $x_1^{j_1}\cdot\ldots\cdot
x_k^{j_k}$. When the monom is evaluated, we get
$$\pi^{\log{r_1^{j_1}}+\ldots+\log{r_k^{j_k}}}=
\pi^{\log{\left(r_1^{j_1}\cdot\ldots\cdot r_k^{j_k}\right)}}$$
Now it is easy to see that if an integer combination of such terms is
zero, then the combination must be either trivial or sum to $0$ and
all monoms must be equal. Both cases imply a polynomial identically
equal to zero. So, any non-trivial integer polynomial evaluated at $x$
must be non-zero.
So, having this result in hand, now it is straightforward to calculate
higher moments of normal distribution. Here we define a container,
which does the job. In its constructor, we simply calculate Kronecker
powers of $v$ and apply $F_n$ to $\otimes^nv$. $F_n$ is, in fact, a
set of all equivalences in sense of class |Equivalence| over $2n$
elements, having $n$ classes each of them having exactly 2 elements. */
#ifndef NORMAL_MOMENTS_H
#define NORMAL_MOMENTS_H
#include "t_container.hh"
class UNormalMoments : public TensorContainer<URSingleTensor>
{
public:
UNormalMoments(int maxdim, const TwoDMatrix &v);
private:
void generateMoments(int maxdim, const TwoDMatrix &v);
static bool selectEquiv(const Equivalence &e);
};
class FNormalMoments : public TensorContainer<FRSingleTensor>
{
public:
FNormalMoments(const UNormalMoments &moms);
};
#endif