This file contains supplementary data for Ivo Alberto Valente Tavares, “Uncertainty quantification with a Gaussian Process Prior : an example from macroeconomics,” http://hdl.handle.net/10400.5/21444.

Ireland 2004

Ivo Tavares

Model Discrepancy

◼

Learning (no mean, no R)

Auxiliary Functions I (vectorize by columns, inverse of vectorize by columns)

In[]:=

(*Thisvectorizeisequaltovec().Forvech()wesimplyuseArrayReshape[].FortheXandYvariablesweusevech.Forallotherweusevec.*)vectorize[x_,len_]:=ArrayReshape[Transpose[x],{len,1}];invVectorize[x_,{r_,c_}]:=Transpose[ArrayReshape[x,{c,r}]];(*Wechangetheorderfrom{r,c}to{c,r},becausetheArrayReshapeisbyrowsinsteadofbycolumns.TheinverseofvechisArrayReshape[,{r,c}].I'llonlyuseArrayReshapeifIwanttorebuildanoriginalmatrixwhichwasvectorizedwithvech*)

Covariance Function

In[]:=

k[ν_,l_,x_,y_]:=Module[{},(*EstafunçãodáproblemasquandocalculadacomNorm[x-y]=0...Temdesefazerolimite,quedá1.*)If[Norm[x-y]0,1,(*2^(1-ν)/Gamma[ν]*(Norm[x-y]*Sqrt[2*ν]/l)^ν*BesselK[ν,Norm[x-y]*Sqrt[2*ν]/l]ThisisforageneralMaternfunction*)Exp[-Norm[x-y]/l]]];(*Kint=K_{d1,d2}(X1,X2)*)Kint[X1_,X2_,l_]:=Module[{},Table[Table[k[1/2,l,X1〚i〛,X2〚j〛],{j,1,Length[X2]}],{i,1,Length[X1]}](*we'reusingtherowsasvectors.Inreality,we'llsaveXintheprogramasthetransposeoftheXinthetext*)];KK[tildeX_,M1_]:=Module[{A},A=Kint[tildeX,tildeX];(*ArrayFlatten[Table[Table[A,{i,1,dim}],{j,1,dim}]]It'sfasterwiththekroneckerproduct*)KroneckerProduct[M1,A]]

RBC Model

In[]:=

rbcACmat[{β_,δ_,θ_,ρ_,aPar_,γ_}]:=Module[{L1,L2,K11,K12,K22,S1,S2,S3,S4,S5,S6,MatPi,W,U,Cmat,CmatStar},η=1.0051;L1=(η-β*(1-δ))/(β*η*θ^2);L2=ρ*(η-β*(1-δ))/(η-β*(1-θ)*(1-δ));K11=(η-β*(1-θ)*(1-δ))/(β*η*θ);K12=(β*η*θ^2-η+β*(1-θ^2)*(1-δ))/(β*η*θ^2);K22=η*θ/(η-β*(1-θ)*(1-δ));S1=(K22-K11)/K12;S2=((K22-K11)*L1-K12*L2)/(K12*(K11-ρ));S3=K22;S4=K12*L2/(K22-K11)+((K22-ρ)*((K22-K11)*L1-K12*L2))/((K22-K11)*(K11-ρ));S5={1-(1-θ)/θ*S1,S1+(η/β-1+δ)/(θ^2*(η-1+δ))*(θ-S1),1-1/θ*S1};S6={1/θ-(1-θ)/θ*S2,S2+(η/β-1+δ)/(θ^2*(η-1+δ))*(1-S2),1/θ-1/θ*S2};U=Transpose[{Flatten[{S5,S1}],Flatten[{S6,S2}]}];MatPi={{S3,S4},{0,ρ}};(*W={{0},{1}};*)ysteady=aPar^(1/(1-θ))*(θ/(η/β-1+δ))^(θ/(1-θ))*(1-θ)/γ*(1-θ*(η-1+δ)/(η/β-1+δ))^-1;csteady=(1-θ*(η-1+δ)/(η/β-1+δ))*ysteady;hstead=(1-θ)/γ*(1-θ*(η-1+δ)/(η/β-1+δ))^-1;Cmat=Drop[U,{2}];CmatStar={{Log[ysteady]},{Log[hstead]},{Log[csteady]}};{MatPi,Cmat,CmatStar}]

Auxiliary Functions II (reparametrizations, sparse matrices, etc)

Some functions will not be used in the final simulations, namely those that build the matrix derivatives. Check the Appendix on the Phd - Thesis. However, I' ve decided to keep them for completeness.

In[]:=

(**)spmatSigma[T_,dim_]:=Module[{},idT=IdentityMatrix[T,SparseArray];idDim=IdentityMatrix[dim,SparseArray];KroneckerProduct[Flatten[Table[KroneckerProduct[idDim,Transpose[{idT[[i]]}]],{i,1,T}],1],idDim]];spmatL[M1_,T_]:=Module[{},KroneckerProduct[IdentityMatrix[T,SparseArray],Flatten[Table[KroneckerProduct[IdentityMatrix[T,SparseArray],Transpose[{M1[[;;,i]]}]],{i,1,Length[M1]}],1]]];funcHM[KM_,dim_,T_]:=Module[{},KroneckerProduct[Flatten[Table[KroneckerProduct[IdentityMatrix[dim,SparseArray],Transpose[{KM[[;;,i]]}]],{i,T}],1],IdentityMatrix[dim,SparseArray]]];logit[u_]:=Module[{x},x=SetPrecision[u,50];Log[x/(1-x)]];(*melhoraraprecisãodocálculonestasfunções*)invLogit[u_]:=Module[{x},x=SetPrecision[u,50];E^x/(1+E^x)];linTanTransf[u_]:=Module[{x},x=SetPrecision[u,50];Tan[x*Pi/2]];invLinTanTransf[u_]:=Module[{x},x=SetPrecision[u,50];ArcTan[x]/(Pi/2)];reparametrizationMatern[L_]:=Log[L];reparametrizationMaternInv[tilL_]:=Exp[tilL];reparametrizationF[{β_,δ_,θ_,ρ_,aPar_,γ_}]:=Module[{},{logit[β],logit[δ],logit[θ],logit[ρ],Log[aPar],Log[γ]}];(*Transformsoriginalthetasintotildethetas*)reparametrizationFInv[{tilβ_,tilδ_,tilθ_,tilρ_,tilaPar_,tilγ_}]:=Module[{},{invLogit[tilβ],invLogit[tilδ],invLogit[tilθ],invLogit[tilρ],Exp[tilaPar],Exp[tilγ]}];(*Transformstildethetasintooriginalthetas*)reparametrizationFQS[{Q_,Sigma_}]:=Module[{},If[Dimensions[Q]=={2,2},{DiagonalMatrix[Log[Diagonal[Q]]],DiagonalMatrix[Log[Diagonal[Sigma]]]},{Log[Q],Log[Sigma]}]];(*TransformsQandSigmaintotildeQandSigma.Workseitherwithmaindiagonalsorthewholematrices.*)reparametrizationFQSInv[{tilQ_,tilSigma_}]:=Module[{},If[Dimensions[tilQ]=={2,2},{DiagonalMatrix[E^Diagonal[tilQ]],DiagonalMatrix[E^Diagonal[tilSigma]]},{E^tilQ,E^tilSigma}]];(*TransformstildeQandSigmaintooriginalQandSigma.Worksbothwithjustthemaindiagonalorthewholematrices*){Astruct,Cmatstruct,CmatStarstruct}=rbcACmat[reparametrizationFInv[{ntilβ,ntilδ,ntilθ,ntilρ,ntilaPar,ntilγ}]];(*A,CmateCmatStarsãoreparametrizadosàcustadostileconomicparameters*)(*UsarFullImplifyparasimplificarfórmulaseseremmaisrápidasacalcular.*)derAstruct[{tilβ_,tilδ_,tilθ_,tilρ_,tilaPar_,tilγ_}]:=Table[D[Astruct[[i,j]],{{ntilβ,ntilδ,ntilθ,ntilρ,ntilaPar,ntilγ}}],{i,1,2},{j,1,2}]/.{ntilβtilβ,ntilδtilδ,ntilθtilθ,ntilρtilρ,ntilaPartilaPar,ntilγtilγ};(*DerivamosAemrelaçãoaosparâmetrostildetheta*)derCmatstruct[{tilβ_,tilδ_,tilθ_,tilρ_,tilaPar_,tilγ_}]:=Table[D[Cmatstruct[[i,j]],{{ntilβ,ntilδ,ntilθ,ntilρ,ntilaPar,ntilγ}}],{i,1,3},{j,1,2}]/.{ntilβtilβ,ntilδtilδ,ntilθtilθ,ntilρtilρ,ntilaPartilaPar,ntilγtilγ};(*DerivamosCmatemrelaçãoaosparâmetrostildetheta*)derCmatStarstruct[{tilβ_,tilδ_,tilθ_,tilρ_,tilaPar_,tilγ_}]:=Table[D[CmatStarstruct[[i,j]],{{ntilβ,ntilδ,ntilθ,ntilρ,ntilaPar,ntilγ}}],{i,1,3},{j,1,1}]/.{ntilβtilβ,ntilδtilδ,ntilθtilθ,ntilρtilρ,ntilaPartilaPar,ntilγtilγ};(*DerivamosCmatStaremrelaçãoaosparâmetrostildetheta*)derKM[X1T_,tilL_]:=Module[{},Table[Table[k[1/2,Exp[tilL],X1T〚i〛,X1T〚j〛]*Norm[X1T〚i〛-X1T〚j〛]*Exp[-tilL],{j,1,Length[X1T]}],{i,1,Length[X1T]}]];upperTriangular[v_,n_]:=Module[{i=0,j=0},Array[If[#>#2,0,v[[++i]]]&,{n,n}]];(*Tocreateanuppertriangularmatrixfromavector.WeshallusethisindrawPropforM1matrix*)

Auxiliary Functions III (Normal density,positivizeMatrix, etc)

In[]:=

(*transitionDrawPGandcondNormalwillbeusedinthePGASMKalgorithm*)transitionDrawPG[xt1_,A_,Q_]:=RandomVariate[MultinormalDistribution[Flatten[A.xt1],Q]];positivizeMatrix[X_]:=Module[{res},res=X;res=0.5*(res+Transpose[res]);(*ThismethodonlyworksforHermitianmatrices*)If[Not[PositiveDefiniteMatrixQ[res]]||MatrixRank[res]≠Length[res],auxsystDP=Eigensystem[res];duDP=DiagonalMatrix[auxsystDP[[1]]];wDP=ReplacePart[duDP,{i_,i_}/;duDP[[i,i]]≤00];lenRes=Length[res];powerDP=30;While[Not[PositiveDefiniteMatrixQ[res]]||MatrixRank[res]≠Length[res],res=Transpose[auxsystDP[[2]]].(wDP+IdentityMatrix[lenRes]*$MinMachineNumber*2^(-powerDP)).auxsystDP[[2]];res=0.5*(res+Transpose[res]);(*MachinePrecisionconsiderationsmaketheabovemandatory!Nottheoreticallyneeded.*)powerDP--;];];res];logMultinormalDens[x_,mean_,var2_]:=Module[{t,det,var},(*Theoretically,I'vemadesurevar2isPDandSym,however,duetomachine-precisionconsiderations,wemustmakesureit'ssymmetric.Otherwisewewouldhaveanerrorwarningfrommathematica*)var=positivizeMatrix[var2];t=Quiet[0.5*(x-mean).LinearSolve[var,(x-mean),Method"Cholesky"]];det=(2*Pi)^Length[var]*Det[var]^(0.5);-t-Log[det]];logCondNormal[vecy1T_,X1T_,t1_,t2_,M1PG_,tilL_,CmatPG_,CmatStarPG_,SigmaPG_]:=Module[{},(*{fvecy1T,vecy1tminus1,vecytT,mean,mean1tminus1,meantT,varMat,varMatTT,varMat1t,varMat1tT,meancond,varcond,fvy}*)(*N(y1|y2)comdim(y1)=3*t1edim(y1)+dim(y2)=3*t2*)(*vecy1Ttemdeserflattened*)fvecy1T=Flatten[vecy1T];vecy1tminus1=fvecy1T[[1;;3*t1]];(*(t1-1)*)vecytT=fvecy1T[[3*t1+1;;3*t2]];mean=Flatten[KroneckerProduct[IdentityMatrix[t2],CmatPG].ArrayReshape[X1T,{2*t2,1}]+KroneckerProduct[ConstantArray[1,t2],IdentityMatrix[3]].CmatStarPG];(*WemustuseArrayReshape,sinceX1TwasmadeusingobservationsasROWS,Notcolumns*)mean1tminus1=mean[[1;;3*t1]];meantT=mean[[3*t1+1;;]];X1t2=X1T[[1;;t2]];varMat=positivizeMatrix[KroneckerProduct[Kint[X1t2,X1t2,Exp[tilL]],M1PG]+KroneckerProduct[IdentityMatrix[t2],SigmaPG]];varMatTT=varMat[[3*t1+1;;,3*t1+1;;]];varMat1t=varMat[[1;;3*t1,1;;3*t1]];varMat1tT=varMat[[1;;3*t1,3*t1+1;;]];meancond=Quiet[Flatten[meantT+Transpose[varMat1tT].LinearSolve[varMat1t,(vecy1tminus1-mean1tminus1)]]];varcond=Quiet[varMatTT-Transpose[varMat1tT].LinearSolve[varMat1t,varMat1tT]];(*Theoreticallythenextlinesofcodeshouldn'tbeneeded,butduetonumericalmachineaccuracy,wesometimesneedthis*)varcond=positivizeMatrix[varcond];t=Quiet[0.5*(vecytT-meancond).LinearSolve[varcond,(vecytT-meancond),Method"Cholesky"]];(*Quietquietsnon-criticalwarnings,likenumericallyunstablematrices.*)det=(2*Pi)^Length[varcond]*Det[varcond]^(0.5);(*Theoreticallytheabovelineofcodeshouldn'tneedAbs,butduetonumericalmachineaccuracy,wesometimesneedthis*)-t-Log[det]];

Mean for Econ

In[]:=

newMeanEcon[vecy1T_,vecx2T_,vecx1T1_,vecx1T_,X1T_,T_,KM_,{M1_,{tilβ_,tilδ_,tilθ_,tilρ_,tilaPar_,tilγ_},tilL_,{tilQMat_,tilSigmaMat_}},s1_]:=Module[{},Tilparam={tilβ,tilδ,tilθ,tilρ,tilaPar,tilγ}+s1*RandomReal[{-1,1},6];Tilparam];

Mean for Vars

In[]:=

newMeanVars[vecy1T_,vecx2T_,vecx1T1_,vecx1T_,X1T_,T_,KM_,{M1_,{tilβ_,tilδ_,tilθ_,tilρ_,tilaPar_,tilγ_},tilL_,{tilQMat_,tilSigmaMat_}},s1_]:=Module[{},resTilVars={tilQMat+s1*DiagonalMatrix[RandomReal[{-1,1},2]],tilSigmaMat+s1*DiagonalMatrix[RandomReal[{-1,1},3]]};resTilVars];

Mean for tilL

In[]:=

newMeanL[vecy1T_,vecx2T_,vecx1T1_,vecx1T_,X1T_,T_,KM_,{M1_,{tilβ_,tilδ_,tilθ_,tilρ_,tilaPar_,tilγ_},tilL_,{tilQMat_,tilSigmaMat_}},s1_]:=Module[{},newtiL=tilL+s1*RandomReal[{-1,1}];newtiL];

Mean for M1(in the proposal we adjust for the d.f.)

In[]:=

newMeanM1[vecy1T_,vecx2T_,vecx1T1_,vecx1T_,X1T_,T_,KM_,{M1_,{tilβ_,tilδ_,tilθ_,tilρ_,tilaPar_,tilγ_},tilL_,{tilQMat_,tilSigmaMat_}},s1_]:=Module[{},CovMat=IdentityMatrix[3];(*RandomVariate[InverseWishartMatrixDistribution[3,IdentityMatrix[3]]];*)newM1=M1+s1*s1M1*RandomVariate[InverseWishartMatrixDistribution[nuM1,CovMat]];newM1];

Marginal Likelihood

In[]:=

logMarginalLikelihood[vecy1T_,X1T_,T_,{M1_,{tilβ_,tilδ_,tilθ_,tilρ_,tilaPar_,tilγ_},tilL_,{tilQMat_,tilSigmaMat_}}]:=Module[{},vecx1TLik=ArrayReshape[X1T,{Times@@Dimensions[X1T],1}];(*WemustuseArrayReshapesinceX1TishasineachROWanobservation*){ALik,CmatLik,CmatStarLik}=rbcACmat[reparametrizationFInv[{tilβ,tilδ,tilθ,tilρ,tilaPar,tilγ}]];{QLik,SigmaLik}=reparametrizationFQSInv[{tilQMat,tilSigmaMat}];If[T==1,KMLik={{1}};resLik=logMultinormalDens[X1T[[1]],{0,0},IdentityMatrix[2]]+logMultinormalDens[Flatten[vecy1T],Flatten[KroneckerProduct[IdentityMatrix[T],CmatLik].Flatten[vecx1TLik]+KroneckerProduct[ConstantArray[1,T],IdentityMatrix[3]].CmatStarLik],KroneckerProduct[KMLik,M1]+KroneckerProduct[IdentityMatrix[T],SigmaLik]]];If[T≥2,x2TLik=Drop[X1T,1];x1T1Lik=Drop[X1T,-1];(*WemustuseArrayReshapesinceX1ThasineachROWanobservation*)vecx2TLik=ArrayReshape[x2TLik,{Times@@Dimensions[x2TLik],1}];vecx1T1Lik=ArrayReshape[x1T1Lik,{Times@@Dimensions[x1T1Lik],1}];KMLik=Kint[X1T,X1T,Exp[tilL]];resLik=logMultinormalDens[X1T[[1]],{0,0},IdentityMatrix[2]]+logMultinormalDens[Flatten[vecx2TLik],Flatten[(KroneckerProduct[IdentityMatrix[T-1],ALik].vecx1T1Lik)],KroneckerProduct[IdentityMatrix[T-1],QLik]]+logMultinormalDens[Flatten[vecy1T],Flatten[KroneckerProduct[IdentityMatrix[T],CmatLik].vecx1TLik+KroneckerProduct[ConstantArray[1,T],IdentityMatrix[3]].CmatStarLik],KroneckerProduct[KMLik,M1]+KroneckerProduct[IdentityMatrix[T],SigmaLik]];];resLik];

Prior(Options for M1: Chang Eaves, Wishart, Inv-Wishart)

In[]:=

(*DefinemultivariategammaforthedensityofInverseWishart*)multiΓ[p_,a_]:=π^(p*(p-1)/4)Product[Gamma[a+(1-j)/2],{j,1,p}]

In[]:=

(*DensityfunctionofinverseWishart*)logDensInvWish[Σ_,ν_,X2_]:=Module[{len,eigen,lbnugget,cond,logpropconst},(*νmustbegreaterthannumberoflines/columnsofΣ*)(*machine-precisionconsiderations*)len=Length[Σ];Xdens=0.5*(X2+Transpose[X2]);cond=0.5*Tr[Σ.Inverse[Xdens]];logpropconst=Log[Abs[Det[Σ]]^(ν/2)]-Log[(2^(ν*len/2)*multiΓ[len,ν/2])]+Log[Norm[Det[Xdens]]^(-(ν+len+1)/2)];(*If[-cond+Log[propconst]<10.Log[$MinMachineNumber],0(*Weneedthis,otherwisewe'llhaveunderflowproblems,comingfromExp[-cond]*),propconst*Exp[-cond]]*)-cond+logpropconst];

In[]:=

(*KerneldensityforChangandEaves'sObjectivePrior*)logChangEaves[X_]:=Module[{len,cond,eigen,lbnugget},len=Length[X];Xdens=0.5*(X+Transpose[X]);Log[Norm[Det[Xdens]]^(-0.5*(len+1))]+Log[Norm[Det[IdentityMatrix[len]+Xdens*Inverse[Xdens]]]^(-0.5)](*ThedensityisPROPORTIONALtotheresult*)];absDetJacob[{tilβ_,tilδ_,tilθ_,tilρ_,tilaPar_,tilγ_},tilL_,{tilQMat_,tilSigmaMat_}]:=Module[{},Abs[E^tilβ/(1+E^tilβ)^2*E^tilδ/(1+E^tilδ)^2*E^tilθ/(1+E^tilθ)^2*E^tilρ/(1+E^tilρ)^2*Exp[tilaPar+tilγ+tilL]*Apply[Times,E^Diagonal[tilQMat]]*Apply[Times,E^Diagonal[tilSigmaMat]]]];(*AbsoluteValueofJacobian'sDeterminant*)logPriorTheta[{M1_,{tilβ_,tilδ_,tilθ_,tilρ_,tilaPar_,tilγ_},tilL_,{tilQMat_,tilSigmaMat_}}]:=Module[{},{βprior,δprior,θprior,ρprior,aParprior,γprior}=reparametrizationFInv[{tilβ,tilδ,tilθ,tilρ,tilaPar,tilγ}];{Qprior,Sigmaprior}=reparametrizationFQSInv[{tilQMat,tilSigmaMat}];(*Ifarhodrawisnegative,thenwe'llgetazerofromtheprior.*)Log[PDF[BetaDistribution[9.9,1.7],βprior]*PDF[BetaDistribution[1,10],δprior]*PDF[BetaDistribution[5,10],θprior]*PDF[BetaDistribution[9.9,1.7],ρprior]*PDF[GammaDistribution[10,0.5],aParprior]*PDF[GammaDistribution[0.25,0.1],γprior]*PDF[GammaDistribution[2,2],Exp[tilL]]*PDF[InverseGammaDistribution[2.5,0.5],Qprior[[1,1]]]*PDF[InverseGammaDistribution[2.5,0.5],Qprior[[2,2]]]*PDF[InverseGammaDistribution[2.5,0.5],Sigmaprior[[1,1]]]*PDF[InverseGammaDistribution[2.5,0.5],Sigmaprior[[2,2]]]*PDF[InverseGammaDistribution[2.5,0.5],Sigmaprior[[3,3]]]*absDetJacob[{tilβ,tilδ,tilθ,tilρ,tilaPar,tilγ},tilL,{tilQMat,tilSigmaMat}]]+logDensInvWish[IdentityMatrix[3],15,M1]];

Draws for Proposals and Pdfs

In[]:=

drawPropEcon[mean_,s2_]:=Module[{draw},scaleMatrixEcon=IdentityMatrix[6];(*If[Select[mean,Element[#,Reals]&]≠mean,mean2={tilβ,tilδ,tilθ,tilρ,tilaPar,tilγ};,mean2=mean;];*)(*scaleMatrixEcon[[4,4]]=5;*)draw=RandomVariate[MultivariateTDistribution[mean,s2*scaleMatrixEcon,3]];draw];pdfPropEcon[x_,mean_,s2_]:=Module[{pdf},pdf=PDF[MultivariateTDistribution[mean,s2*scaleMatrixEcon,3],x];(*Check[pdf=PDF[MultivariateTDistribution[mean,s2*scaleMatrixEcon,3],x],Print["x = ",InputForm[x]];Print["mean = ",InputForm[mean]];pdf=PDF[MultivariateTDistribution[mean,s2*scaleMatrixEcon,3],x]];*)pdf];

In[]:=

drawPropVars[{meanQ_,meanSigma_},s2_]:=Module[{},scaleMatrixVars=IdentityMatrix[5];mean=Flatten[{Diagonal[meanQ],Diagonal[meanSigma]}];tSampleDP=RandomVariate[MultivariateTDistribution[mean,s2*scaleMatrixVars,3]];newtilQDP=DiagonalMatrix[tSampleDP[[1;;2]]];newtilSigmaDP=DiagonalMatrix[tSampleDP[[3;;5]]];{newtilQDP,newtilSigmaDP}(*drawPropoutputsthetildeparametersfortheeconomicstructuralparameters,andalsotildeQandtildeSigma.Everythingasmatrices,exceptresnewstrucwhichgoesoutasavector,andresnewLgpasanumber*)];pdfPropVars[{xQ_,xSigma_},{meanQ_,meanSigma_},s2_]:=Module[{pdf},pdf=PDF[MultivariateTDistribution[Flatten[{Diagonal[meanQ],Diagonal[meanSigma]}],s2*scaleMatrixVars,3],Flatten[{Diagonal[xQ],Diagonal[xSigma]}]];pdf];

In[]:=

drawPropL[mean_,s2_]:=Module[{},RandomVariate[StudentTDistribution[mean,s2,3]]];pdfPropL[x_,mean_,s2_]:=Module[{},PDF[StudentTDistribution[mean,s2,3],x]];

In[]:=

drawPropM1[mean_,s2_]:=Module[{},(*nu=5;*)matrixM1Prop=(*mean;*)positivizeMatrix[mean*(nuM1+Length[mean]+1)];(*Weneedthis,sincesometimesmultiplyingbynu+Length[mean]+1willreturnanon-SPDmatrixduetomachineprecisionconsiderations...*)RandomVariate[InverseWishartMatrixDistribution[nuM1,matrixM1Prop]]];logPdfPropM1[x_,mean_,s2_]:=Module[{},(*nu=5;*)logDensInvWish[mean*(nuM1+Length[mean]+1),nuM1,x]]

Auxiliary Functions 4 (Acceptance Step for MH Steps)

In[]:=

AcceptanceStep[vecy1T_,X1T_,T_,thetaMinusList_,thetaStarList_,newMeanGivenThetaMinus1_,newMeanGivenThetaStarList_,s1_,s2_,option_]:=Module[{},(*thetaMinusList={{tilβ,tilδ,tilθ,tilρ,tilaPar,tilγ},tilL,tilQMat};*)(*option1,EconStepoption2,LStepoption3,Q,SigmaStepoption4,M1step*)logdenominatorprior=logPriorTheta[thetaMinusList];logdenominatorproposal=Which[option1,Log[pdfPropEcon[thetaStarList[[2]],newMeanGivenThetaMinus1,s2]],option2,Log[pdfPropL[thetaStarList[[3]],newMeanGivenThetaMinus1,s2]],option3,Log[pdfPropVars[thetaStarList[[4]],newMeanGivenThetaMinus1,s2]],option4,logPdfPropM1[thetaStarList[[1]],newMeanGivenThetaMinus1,s2]];(*thetaStarList={M1,drawPropEcon[newMeanEconGivenThetaMinus1,s2],tilL,{tilQMat,tilSigmaMat}};*)logdenominatorlikelihood=logMarginalLikelihood[vecy1T,X1T,T,thetaMinusList];logdenominator=logdenominatorprior+logdenominatorproposal+logdenominatorlikelihood;(*---Numerator---*)lognumeratorprior=logPriorTheta[thetaStarList];lognumeratorlik=logMarginalLikelihood[vecy1T,X1T,T,thetaStarList];lognumeratorproposal=Which[option1,Log[pdfPropEcon[thetaMinusList[[2]],newMeanGivenThetaStarList,s2]],option2,Log[pdfPropL[thetaMinusList[[3]],newMeanGivenThetaStarList,s2]],option3,Log[pdfPropVars[thetaMinusList[[4]],newMeanGivenThetaStarList,s2]],option4,logPdfPropM1[thetaMinusList[[1]],newMeanGivenThetaStarList,s2]];(*thetaStarList={M1,drawPropEcon[newMeanEconGivenThetaMinus1,s2],tilL,{tilQMat,tilSigmaMat}};*)logAcceptprob=Min[{0,(lognumeratorprior+lognumeratorproposal+lognumeratorlik-logdenominator)}];(*Log[1]=0*)u=RandomVariate[UniformDistribution[{0,1}],1][[1]];If[Log[u]≤logAcceptprob,res=thetaStarList;counterMH=1;,counterMH=0;res=thetaMinusList;];{res,counterMH}];

M-H step Econ

In[]:=

MHstepEcon[vecy1T_,X1T_,T_,dim_,{M1_,{tilβ_,tilδ_,tilθ_,tilρ_,tilaPar_,tilγ_},tilL_,{tilQMat_,tilSigmaMat_}},s1_,s2_]:=Module[{},(*Print[" - - Entrando no MH - - "];*)(*x2TMH=Drop[X1T,1];x1T1MH=Drop[X1T,-1];(*WeuseArrayReshapesincethematrixX1TismadeofROWSofobservations*)vecx2TMH=ArrayReshape[x2TMH,{Times@@Dimensions[x2TMH],1}];vecx1TMH=ArrayReshape[X1T,{Times@@Dimensions[X1T],1}];vecx1T1MH=ArrayReshape[x1T1MH,{Times@@Dimensions[x1T1MH],1}];*)(*TheVectorizationsabovewillbeintheMHGlobaltobemoreefficient.*)KMMH=Kint[X1T,X1T,Exp[tilL]];{AMH,CmatMH,CmatStarMH}=rbcACmat[reparametrizationFInv[{tilβ,tilδ,tilθ,tilρ,tilaPar,tilγ}]];{QMH,SigmaMH}=reparametrizationFQSInv[{tilQMat,tilSigmaMat}];thetaMinusList={M1,{tilβ,tilδ,tilθ,tilρ,tilaPar,tilγ},tilL,{tilQMat,tilSigmaMat}};newMeanEconGivenThetaMinus1=newMeanEcon[vecy1T,vecx2TMH,vecx1T1MH,vecx1TMH,X1T,T,KMMH,{M1,{tilβ,tilδ,tilθ,tilρ,tilaPar,tilγ},tilL,{tilQMat,tilSigmaMat}},s1];(*meanforq(θ_*|θ[n-1],X1:T,y1:T)*)thetaStarList={M1,drawPropEcon[newMeanEconGivenThetaMinus1,s2],tilL,{tilQMat,tilSigmaMat}};(*thetastarlististhecandidatethetadrawasalistoflists.We'llneedtotransformitintoapropervectorusingonlytheuppertriangularmatrixofM1,themaindiagonalofQandSigma*)KMMHThetaStar=KMMH;(*InthisMHstep,KMMH=Kint[X1T,X1T,Exp[thetaStarList〚3〛]];*)newMeanEconGivenThetaStarList=newMeanEcon[vecy1T,vecx2TMH,vecx1T1MH,vecx1TMH,X1T,T,KMMHThetaStar,thetaStarList,s1];(*meanforq(θ[n-1]|θ_*,X1:T,y1:T)*)(*thetaStaristhetheta_*vectorwhichwillenterintheproposaldensity*)AcceptanceStep[vecy1T,X1T,T,thetaMinusList,thetaStarList,newMeanEconGivenThetaMinus1,newMeanEconGivenThetaStarList,s1,s2,1](*AcceptanceStep[vecy1T_,X1T_,T_,thetaMinusList_,thetaStarList_,newMeanGivenThetaMinus1_,newMeanGivenThetaStarList_,s1_,s2_,option_]*)];

M-H Step L

In[]:=

MHstepL[vecy1T_,X1T_,T_,dim_,{M1_,{tilβ_,tilδ_,tilθ_,tilρ_,tilaPar_,tilγ_},tilL_,{tilQMat_,tilSigmaMat_}},s1_,s2_]:=Module[{},KMMH=Kint[X1T,X1T,Exp[tilL]];{AMH,CmatMH,CmatStarMH}=rbcACmat[reparametrizationFInv[{tilβ,tilδ,tilθ,tilρ,tilaPar,tilγ}]];{QMH,SigmaMH}=reparametrizationFQSInv[{tilQMat,tilSigmaMat}];thetaMinusList={M1,{tilβ,tilδ,tilθ,tilρ,tilaPar,tilγ},tilL,{tilQMat,tilSigmaMat}};newMeanLGivenThetaMinus1=newMeanL[vecy1T,vecx2TMH,vecx1T1MH,vecx1TMH,X1T,T,KMMH,{M1,{tilβ,tilδ,tilθ,tilρ,tilaPar,tilγ},tilL,{tilQMat,tilSigmaMat}},s1];(*meanforq(θ_*|θ[n-1],X1:T,y1:T)*)thetaStarList={M1,{tilβ,tilδ,tilθ,tilρ,tilaPar,tilγ},drawPropL[newMeanLGivenThetaMinus1,s2],{tilQMat,tilSigmaMat}};(*thetastarlististhecandidatethetadrawasalistoflists.We'llneedtotransformitintoapropervectorusingonlytheuppertriangularmatrixofM1,themaindiagonalofQandSigma*)KMMHThetaStar=Kint[X1T,X1T,Exp[thetaStarList〚3〛]];(*InthisMHstepKMMHisdifferent,becausethetaStarList〚3〛isdifferentfromtilL*)newMeanLGivenThetaStarList=newMeanL[vecy1T,vecx2TMH,vecx1T1MH,vecx1TMH,X1T,T,KMMHThetaStar,thetaStarList,s1];(*meanforq(θ[n-1]|θ_*,X1:T,y1:T)*)AcceptanceStep[vecy1T,X1T,T,thetaMinusList,thetaStarList,newMeanLGivenThetaMinus1,newMeanLGivenThetaStarList,s1,s2,2]];

M-H Step Vars

In[]:=

MHstepVars[vecy1T_,X1T_,T_,dim_,{M1_,{tilβ_,tilδ_,tilθ_,tilρ_,tilaPar_,tilγ_},tilL_,{tilQMat_,tilSigmaMat_}},s1_,s2_]:=Module[{},KMMH=Kint[X1T,X1T,Exp[tilL]];{AMH,CmatMH,CmatStarMH}=rbcACmat[reparametrizationFInv[{tilβ,tilδ,tilθ,tilρ,tilaPar,tilγ}]];{QMH,SigmaMH}=reparametrizationFQSInv[{tilQMat,tilSigmaMat}];thetaMinusList={M1,{tilβ,tilδ,tilθ,tilρ,tilaPar,tilγ},tilL,{tilQMat,tilSigmaMat}};newMeanVarsGivenThetaMinus1=newMeanVars[vecy1T,vecx2TMH,vecx1T1MH,vecx1TMH,X1T,T,KMMH,{M1,{tilβ,tilδ,tilθ,tilρ,tilaPar,tilγ},tilL,{tilQMat,tilSigmaMat}},s1];(*meanforq(θ_*|θ[n-1],X1:T,y1:T)*)thetaStarList={M1,{tilβ,tilδ,tilθ,tilρ,tilaPar,tilγ},tilL,drawPropVars[newMeanVarsGivenThetaMinus1,s2]};(*thetastarlististhecandidatethetadrawasalistoflists.We'llneedtotransformitintoapropervectorusingonlytheuppertriangularmatrixofM1,themaindiagonalofQandSigma*)(*Print["2"];*)KMMHThetaStar=KMMH;(*InthisMHstep,KMMH=Kint[X1T,X1T,Exp[thetaStarList〚3〛]];*)newMeanVarsGivenThetaStarList=newMeanVars[vecy1T,vecx2TMH,vecx1T1MH,vecx1TMH,X1T,T,KMMHThetaStar,thetaStarList,s1];(*meanforq(θ[n-1]|θ_*,X1:T,y1:T)*)(*Print["3"];*)AcceptanceStep[vecy1T,X1T,T,thetaMinusList,thetaStarList,newMeanVarsGivenThetaMinus1,newMeanVarsGivenThetaStarList,s1,s2,3]];

M-H Step M1

In[]:=

MHstepM1[vecy1T_,X1T_,T_,dim_,{M1_,{tilβ_,tilδ_,tilθ_,tilρ_,tilaPar_,tilγ_},tilL_,{tilQMat_,tilSigmaMat_}},s1_,s2_]:=Module[{},KMMH=Kint[X1T,X1T,Exp[tilL]];{AMH,CmatMH,CmatStarMH}=rbcACmat[reparametrizationFInv[{tilβ,tilδ,tilθ,tilρ,tilaPar,tilγ}]];{QMH,SigmaMH}=reparametrizationFQSInv[{tilQMat,tilSigmaMat}];thetaMinusList={M1,{tilβ,tilδ,tilθ,tilρ,tilaPar,tilγ},tilL,{tilQMat,tilSigmaMat}};newMeanM1GivenThetaMinus1=newMeanM1[vecy1T,vecx2TMH,vecx1T1MH,vecx1TMH,X1T,T,KMMH,{M1,{tilβ,tilδ,tilθ,tilρ,tilaPar,tilγ},tilL,{tilQMat,tilSigmaMat}},s1];(*meanforq(θ_*|θ[n-1],X1:T,y1:T)*)(*Print[1];*)thetaStarList={drawPropM1[newMeanM1GivenThetaMinus1,s2],{tilβ,tilδ,tilθ,tilρ,tilaPar,tilγ},tilL,{tilQMat,tilSigmaMat}};(*thetastarlististhecandidatethetadrawasalistoflists.We'llneedtotransformitintoapropervectorusingonlytheuppertriangularmatrixofM1,themaindiagonalofQandSigma*)(*Print[2];*)KMMHThetaStar=KMMH;(*InthisMHstepKMMH=Kint[X1T,X1T,Exp[thetaStarList〚3〛]]*)newMeanM1GivenThetaStarList=newMeanM1[vecy1T,vecx2TMH,vecx1T1MH,vecx1TMH,X1T,T,KMMHThetaStar,thetaStarList,s1];(*meanforq(θ[n-1]|θ_*,X1:T,y1:T)*)(*Print["Before acceptstep"];*)AcceptanceStep[vecy1T,X1T,T,thetaMinusList,thetaStarList,newMeanM1GivenThetaMinus1,newMeanM1GivenThetaStarList,s1,s2,4]];

M-H Global

In[]:=

MHstepGlobal[vecy1T_,X1T_,T_,dim_,{M1_,{tilβ_,tilδ_,tilθ_,tilρ_,tilaPar_,tilγ_},tilL_,{tilQMat_,tilSigmaMat_}},s1_,s2_]:=Module[{},x2TMH=Drop[X1T,1];x1T1MH=Drop[X1T,-1];(*WeuseArrayReshapesincethematrixX1TismadeofROWSofobservations*)vecx2TMH=ArrayReshape[x2TMH,{Times@@Dimensions[x2TMH],1}];vecx1TMH=ArrayReshape[X1T,{Times@@Dimensions[X1T],1}];vecx1T1MH=ArrayReshape[x1T1MH,{Times@@Dimensions[x1T1MH],1}];(*ThevectorizationsabovewillbeusedbyalltheMHStepsdefinedabove.*)(*Print["MHstepEcon"];*){resGlobal,acceptEcon}=MHstepEcon[vecy1T,X1T,T,dim,{M1,{tilβ,tilδ,tilθ,tilρ,tilaPar,tilγ},tilL,{tilQMat,tilSigmaMat}},s1,s2];(*Print["MHstepVars"];*){resGlobal,acceptVars}=MHstepVars[vecy1T,X1T,T,dim,resGlobal,s1,s2];(*Print["MHstepL"];*){resGlobal,acceptL}=MHstepL[vecy1T,X1T,T,dim,resGlobal,s1,s2];(*Print["MHstepM1"];*){resGlobal,acceptM1}=MHstepM1[vecy1T,X1T,T,dim,resGlobal,s1,s2];{resGlobal,{acceptEcon,acceptVars,acceptL,acceptM1}}];

PGASMK (Parallelized)

In[]:=

Bayesian Learning and Filtering/Smoothing of SSMs

In[]:=

bayesianLS[time_,vecy1T_,{M1_,{tilβ_,tilδ_,tilθ_,tilρ_,tilaPar_,tilγ_},tilL_,{tilQMat_,tilSigmaMat_}},X1Tinitial_,sampleSize_,numParticles_,stepSize_,varSize_]:=Module[{},(*scaleMatrix=0.5*IdentityMatrix[18];*)outputDim=3;stateSamples=ConstantArray[0,sampleSize];stateSamples[[1]]=X1Tinitial;parameterSamples=ConstantArray[0,sampleSize];parameterSamples[[1]]={M1,{tilβ,tilδ,tilθ,tilρ,tilaPar,tilγ},tilL,{tilQMat,tilSigmaMat}};acceptanceRecord=ConstantArray[0,sampleSize];tildeWRecord=ConstantArray[0,sampleSize];normWRecord=ConstantArray[0,sampleSize];(*thetaStarRecord=ConstantArray[0,sampleSize];*)counterBLS=2;While[counterBLS<=sampleSize,(*Print["Entrando no PGASMK"];*){stateSamples[[counterBLS]],tildeWRecord[[counterBLS-1]],normWRecord[[counterBLS-1]]}=PGASMK[vecy1T,stateSamples[[counterBLS-1]],numParticles,time,outputDim,parameterSamples[[counterBLS-1]]];(*Print["Entrando no MH-Step"];*){parameterSamples[[counterBLS]],acceptanceRecord[[counterBLS-1]]}=MHstepGlobal[vecy1T,stateSamples[[counterBLS]],time,outputDim,parameterSamples[[counterBLS-1]],stepSize,varSize];(*thetaStarRecord[[counterBLS-1]]=thetaStarList;*)(*Print["counterBLS=",counterBLS];*)counterBLS++;];(*pathnameisdefinedbelowatRunCode*)fileDataStates=Compress[stateSamples];Export[pathname<>ToString[time]<>" - EstimationAllSampleStates.m",fileDataStates];fileDataParameters=Compress[parameterSamples];Export[pathname<>ToString[time]<>" - EstimationAllSampleParameters.m",fileDataParameters];fileDataAcceptCounterBLS=Compress[{acceptanceRecord,counterBLS}];Export[pathname<>ToString[time]<>" - EstimationAllSampleAcceptCounterBLS.m",fileDataAcceptCounterBLS];fileDataTildeW=Compress[tildeWRecord];Export[pathname<>ToString[time]<>" - EstimationAllSampleStatesTildeW.m",fileDataTildeW];fileDataNormW=Compress[normWRecord];Export[pathname<>ToString[time]<>" - EstimationAllSampleStatesNormW.m",fileDataNormW];(*fileDataThetaStarRecord=Compress[thetaStarRecord];Export[pathname<>ToString[time]<>" - EstimationAllSampleThetaStarRecord.m",fileDataThetaStarRecord];*){stateSamples,parameterSamples}];

◼

Data(Unfiltered)

Importing Data from File

In[]:=

(*Import["method data.xlsx"];*)

In[]:=

output=Drop[data[[1,All,8]],1];consumption=Drop[data[[1,All,9]],1];hours=Drop[data[[1,All,10]],1];invest=output-consumption;

Data Plots

In[]:=

ListPlot[Drop[Log[output]-Table[i,{i,1,Length[output]}]*Log[1.0051],{1,218-50}]]

Out[]=

In[]:=

ListPlot[Drop[Log[consumption]-Table[i,{i,1,Length[output]}]*Log[1.0051],{1,218-50}]]

Out[]=

Data Treatment and Selection - End Time, Start Time. (See RBC section )

In[]:=

(*Todefinecorrectlythis,wemustobeytheorderintheRBCsectionY_t=(output,hours,consumption)*)intercept=Table[i,{i,1,Length[output]}]*Log[1.0051];{loutput,lhours,lconsumption}={Log[output]-intercept,Log[hours],Log[consumption]-intercept};(*Y1Ttest=Transpose[{loutput,lhours,lconsumption}];endTime=50;startTime=11;loutput=Drop[loutput,{1,Length[output]-endTime}];lconsumption=Drop[lconsumption,{1,Length[output]-endTime}];lhours=Drop[lhours,{1,Length[output]-endTime}];*)Y1T=Transpose[{loutput,lhours,lconsumption}];(*Y1T:eachobservationisstoredinarow.*)

◼

Prediction Functions

Initial Values

In[]:=

initValues[time_]:=Module[{X1,counter,M1,β,δ,θ,ρ,aPar,γ,A,Cmat,CmatStar,Sigma,Q,Lgp},M1={{6.528,1.065,6.565},{1.065,1.607,0.968},{6.565,0.968,6.621}};{β,δ,θ,ρ,aPar,γ}={0.7,0.3,0.6,0.7,2,0.6};{A,Cmat,CmatStar}=rbcACmat[{β,δ,θ,ρ,aPar,γ}];Sigma=DiagonalMatrix[{0.5,0.2,0.3}];Q=DiagonalMatrix[{1.5,0.6}];Lgp=20;X1=ConstantArray[1,time];X1[[1]]=RandomVariate[MultinormalDistribution[{0,0},IdentityMatrix[2]]];counter=2;While[counter≤time,X1[[counter]]=A.X1[[counter-1]]+RandomVariate[MultinormalDistribution[{0,0},Q]];counter++;];{X1,{M1,reparametrizationF[{β,δ,θ,ρ,aPar,γ}],reparametrizationMatern[Lgp],reparametrizationFQS[{Q,Sigma}]}}];

Prediction 1 Step

In[]:=

prediction1Step[stateSamples_,parameterSamples_,time_,vecy1T_,sizeApproxY_,sizeBurnIn_]:=Module[{len},len=Length[stateSamples];(*ForaBurn-In,thepredictionwillonlyusedataaftersizeBurnIn,i.e.sizeBurnIn+1*)econParSampleList=ParallelTable[reparametrizationFInv[parameterSamples[[i,2]]],{i,sizeBurnIn+1,len}];listACCstar=ParallelTable[rbcACmat[econParSampleList[[i]]],{i,1,len-sizeBurnIn}];listQSigma=ParallelTable[reparametrizationFQSInv[parameterSamples[[i,4]]],{i,sizeBurnIn+1,len}];predictedState=Total[ParallelTable[listACCstar[[i-sizeBurnIn,1]].stateSamples[[i,-1]],{i,sizeBurnIn+1,len}]]/(len-sizeBurnIn);(*Thiscomputesx_*.ComputingthemeanstateatT+1,wemusttakecaretouseonlythedataafterBurn-In.*)(*Print["1"];*)KStarX1TList=ParallelTable[KroneckerProduct[Kint[{predictedState},stateSamples[[i]],Exp[parameterSamples[[i,3]]]],parameterSamples[[i,1]]],{i,sizeBurnIn+1,len}];(*KStarX1TListisalistwithK_M(x_*,x_{1:T}[i])\otimesM1[i]ascomponents*)(*Print[2];*)KStarStarList=ParallelTable[KroneckerProduct[Kint[{predictedState},{predictedState},Exp[parameterSamples[[i,3]]]],parameterSamples[[i,1]]],{i,sizeBurnIn+1,len}];(*KStarStarListisalistwithK_M[x_*,x_*]\otimesM1[i]ascomponents*)(*Print[3];*)KX1TList=ParallelTable[KroneckerProduct[Kint[stateSamples[[i]],stateSamples[[i]],Exp[parameterSamples[[i,3]]]],parameterSamples[[i,1]]]+KroneckerProduct[IdentityMatrix[time],listQSigma[[i-sizeBurnIn,2]]],{i,sizeBurnIn+1,len}];(*KX1TListisalistwithK_M[x_{1:T}[i],x_{1:T}[i]]\otimesM1[i]+I_T\otimesSigma[i]ascomponents*)(*Print[4];*)vecX1TList=ParallelTable[ArrayReshape[stateSamples[[i]],{Times@@Dimensions[stateSamples[[i]]],1}],{i,sizeBurnIn+1,len}];meanDifList=ParallelTable[(vecy1T-KroneckerProduct[IdentityMatrix[time],listACCstar[[i,2]]].vecX1TList[[i]]-KroneckerProduct[ConstantArray[1,time],IdentityMatrix[3]].listACCstar[[i,3]]),{i,1,len-sizeBurnIn}];(*meanDifListisalistwithy_{1:T}-(I_T\otimesC)x_{1:T}-C_{*1:T}ascomponents.AllthelistsusedtocreatemeanDifListhavelengthlen-sizeBurnIn*)meanbStarList=ParallelTable[listACCstar[[i,3]]+listACCstar[[i,2]].predictedState+KStarX1TList[[i]].LinearSolve[KX1TList[[i]],meanDifList[[i]]],{i,1,len-sizeBurnIn}];(*meanbStarListgivesalistwith\hat\tildeb_*[i]ascomponents.AllthelistsusedtocreatemeanDifListhavelengthlen-sizeBurnIn*)varbStarList=ParallelTable[positivizeMatrix[KStarStarList[[i]]-KStarX1TList[[i]].LinearSolve[KX1TList[[i]],Transpose[KStarX1TList〚i〛]]],{i,1,len-sizeBurnIn}];(*varbStarListisalistwith\tildeSigma_*[i]ascomponents.Theoretically,thereshouldbenoneedforusingpositivizeMatrixhere,however,duetomachineprecisionconsiderations(roundingerrors),wemustensurethatit'ssymmetricandPD.*)bCatList=Ceiling[RandomVariate[UniformDistribution[{0,1}],len-sizeBurnIn]*(len-sizeBurnIn)];(*Distributionof\tildeb_*isamixtureofdistributions,soweneedto'randomize',withequalprob(hencetheUniformDist),thechoosingofeachcomponentdistribution.*)bStarList=ParallelTable[RandomVariate[MultinormalDistribution[Flatten[meanbStarList〚bCatList〚i〛〛],varbStarList〚bCatList〚i〛〛]],{i,1,len-sizeBurnIn}];(*bStarListisalistof\tildeb_*[i]ascomponents,whichwillbeusedfortheapproximationoftheposteriorpredictivedistribution.*)yCatList=Ceiling[RandomVariate[UniformDistribution[{0,1}],sizeApproxY]*(len-sizeBurnIn)];(*Theposteriorpredictivedistributionofy_*isalsoamixtureofdistributions,soweneedto'randomize',withequalprob(hencetheUniformDist),thechoosingofeachcomponentdistribution.*)yStarList=ParallelTable[RandomVariate[MultinormalDistribution[Flatten[bStarList〚yCatList〚i〛〛],listQSigma〚yCatList〚i〛,2〛]],{i,1,sizeApproxY}];(*prediction1StepreturnsasampledrawnfromanapproximationtotheposteriordistributionofyStar*)fileDataYStarList=Compress[yStarList];(*pathnameisdefinedbelowatRunCode*)Export[pathname<>ToString[time+1]<>" - yStarList.m",fileDataYStarList];];

Prediction N Steps

In[]:=

predictionNSteps[Y1T_,startTime_,endTime_,sizeApproxY_,sampleSize_,numParticles_,stepSize_,varSize_,sizeBurnIn_]:=Module[{},Y1TEnd=Drop[Y1T,{1,Length[Y1T]-endTime}];(*WedropthefirstLength[Y1T]-endTimeobs.*)currentTime=startTime;While[currentTime<endTime,(*wedon'thaveendTime+1observations.*)startValues=initValues[currentTime];(*Print["Y1T treatment"];*)Y1TCurrent=Drop[Y1TEnd,-(endTime-currentTime)];(*WedropthelastendTime-currentTimeobs*)vecy1TCurrent=ArrayReshape[Y1TCurrent,{Times@@Dimensions[Y1TCurrent],1}];(*Print["Y1TCurrent Dimensions = ",Dimensions[Y1TCurrent]];Print["vecy1TCurrent Dimensions = ",Dimensions[vecy1TCurrent]];*)(*Print["BayesianLS"];*){stateSamplesTime,parameterSamplesTime}=bayesianLS[currentTime,vecy1TCurrent,startValues[[2]],startValues[[1]],sampleSize,numParticles,stepSize,varSize];(*Print["prediction1Step"];*)prediction1Step[stateSamplesTime,parameterSamplesTime,currentTime,vecy1TCurrent,sizeApproxY,sizeBurnIn];currentTime++;];];

◼

Run Code for Predictions

In my simulations, I used:
startTime=10;
endTime=15;
sizeApproxY=50000;
sampleSize=15000;
numParticles=20;
stepSize=0.0001;
nuM1=15;
s1M1=1;
varSize=0.1;
sizeBurnIn=5000;

Change the below to your needs, and run the simulation.