(* Content-type: application/vnd.wolfram.mathematica *)

(*** Wolfram Notebook File ***)
(* http://www.wolfram.com/nb *)

(* CreatedBy='Mathematica 10.0' *)

(*CacheID: 234*)
(* Internal cache information:
NotebookFileLineBreakTest
NotebookFileLineBreakTest
NotebookDataPosition[       158,          7]
NotebookDataLength[     25894,        670]
NotebookOptionsPosition[     25272,        644]
NotebookOutlinePosition[     25620,        659]
CellTagsIndexPosition[     25577,        656]
WindowFrame->Normal*)

(* Beginning of Notebook Content *)
Notebook[{
Cell[BoxData[
 RowBox[{
  RowBox[{"Clear", "[", "\"\<Global`*\>\"", "]"}], ";"}]], "Input"],

Cell[CellGroupData[{

Cell[BoxData[{
 RowBox[{
  RowBox[{"k", "=", "3"}], ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"t", " ", "=", " ", "2"}], ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"s1", " ", "=", " ", 
   RowBox[{
    UnderoverscriptBox["\[Sum]", 
     RowBox[{"i", "=", "1"}], "k"], 
    RowBox[{
     RowBox[{"(", 
      RowBox[{"y", "[", "i", "]"}], ")"}], "/", "k"}]}]}], 
  ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{
   RowBox[{"s2", " ", "=", " ", 
    RowBox[{
     UnderoverscriptBox["\[Sum]", 
      RowBox[{"i", "=", "1"}], "k"], 
     RowBox[{
      RowBox[{
       RowBox[{"(", 
        RowBox[{
         RowBox[{"y", "[", "i", "]"}], "-", "s1"}], ")"}], "^", "2"}], "/", 
      "k"}]}]}], ";"}], "\[IndentingNewLine]", 
  RowBox[{"(*", " ", 
   RowBox[{"2", "x2", " ", "matrix"}], " ", "*)"}]}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{
   RowBox[{"x", "[", "1", "]"}], " ", "=", " ", 
   RowBox[{"(", GridBox[{
      {"x111", "x112"},
      {"x121", "x122"}
     }], ")"}]}], ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{
   RowBox[{"x", "[", "2", "]"}], " ", "=", " ", 
   RowBox[{"(", GridBox[{
      {"x211", "x212"},
      {"x221", "x222"}
     }], ")"}]}], ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{
   RowBox[{"x", "[", "3", "]"}], " ", "=", " ", 
   RowBox[{"(", GridBox[{
      {"x311", "x312"},
      {"x321", "x322"}
     }], ")"}]}], ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"m", "=", 
   RowBox[{"(", GridBox[{
      {"m1", "m2"},
      {"m3", "m4"}
     }], ")"}]}], ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{
   RowBox[{"id", " ", "=", " ", 
    RowBox[{"(", GridBox[{
       {"1", "0"},
       {"0", "1"}
      }], ")"}]}], ";"}], "\[IndentingNewLine]", 
  RowBox[{"(*", " ", 
   RowBox[{"implementation", " ", "of", " ", "new", " ", "formula"}], " ", 
   "*)"}]}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"startTime", "=", 
   RowBox[{"AbsoluteTime", "[", "]"}]}], ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"P", "=", "0"}], ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{
   RowBox[{"terms", " ", "=", " ", 
    RowBox[{"MonomialList", "[", 
     RowBox[{"Sum", "[", 
      RowBox[{
       RowBox[{
        FractionBox[
         RowBox[{
          SuperscriptBox[
           RowBox[{"(", 
            RowBox[{"-", "1"}], ")"}], "i"], " ", 
          RowBox[{
           RowBox[{"(", 
            RowBox[{
             RowBox[{"2", "i"}], "-", "1"}], ")"}], "!!"}], 
          RowBox[{"Binomial", "[", 
           RowBox[{"t", ",", 
            RowBox[{"2", "i"}]}], "]"}]}], 
         RowBox[{"Product", "[", 
          RowBox[{
           RowBox[{"k", "+", 
            RowBox[{"2", "j"}], "-", "3"}], ",", 
           RowBox[{"{", 
            RowBox[{"j", ",", "1", ",", "i"}], "}"}]}], "]"}]], 
        SuperscriptBox["s1", 
         RowBox[{"t", "-", 
          RowBox[{"2", "i"}]}]], 
        SuperscriptBox["s2", "i"]}], ",", 
       RowBox[{"{", 
        RowBox[{"i", ",", "0", ",", 
         RowBox[{"Floor", "[", 
          RowBox[{"t", "/", "2"}], "]"}]}], "}"}]}], "]"}], "]"}]}], ";", 
   "  ", 
   RowBox[{"(*", " ", 
    RowBox[{
    "list", " ", "of", " ", "monomial", " ", "terms", " ", "in", " ", "P"}], 
    " ", "*)"}], "\[IndentingNewLine]", "\[IndentingNewLine]", 
   RowBox[{"For", "[", 
    RowBox[{
     RowBox[{"i", "=", "1"}], ",", 
     RowBox[{"i", "\[LessEqual]", 
      RowBox[{"Length", "[", "terms", "]"}]}], ",", 
     RowBox[{"i", "++"}], ",", "\[IndentingNewLine]", 
     RowBox[{
      RowBox[{"term", " ", "=", " ", 
       RowBox[{"terms", "[", 
        RowBox[{"[", "i", "]"}], "]"}]}], ";", "\[IndentingNewLine]", 
      RowBox[{"x1exp", " ", "=", " ", 
       RowBox[{"Exponent", "[", 
        RowBox[{"term", ",", 
         RowBox[{"y", "[", "1", "]"}]}], "]"}]}], ";", "\[IndentingNewLine]", 
      RowBox[{"x2exp", " ", "=", " ", 
       RowBox[{"Exponent", "[", 
        RowBox[{"term", ",", 
         RowBox[{"y", "[", "2", "]"}]}], "]"}]}], ";", "\[IndentingNewLine]", 
      RowBox[{"x3exp", " ", "=", " ", 
       RowBox[{"Exponent", "[", 
        RowBox[{"term", ",", 
         RowBox[{"y", "[", "3", "]"}]}], "]"}]}], ";", "\[IndentingNewLine]", 
      "\[IndentingNewLine]", 
      RowBox[{"exponents", "=", 
       RowBox[{"{", 
        RowBox[{"x1exp", ",", "x2exp", ",", "x3exp"}], "}"}]}], ";", " ", 
      RowBox[{"(*", " ", 
       RowBox[{
       "list", " ", "of", " ", "the", " ", "exponents", " ", "appearing", " ",
         "in", " ", "this", " ", "type", " ", "of", " ", "term"}], " ", 
       "*)"}], "\[IndentingNewLine]", "\[IndentingNewLine]", 
      RowBox[{"type", " ", "=", " ", 
       RowBox[{
        RowBox[{
         RowBox[{"y", "[", "1", "]"}], "^", "x1exp"}], " ", 
        RowBox[{
         RowBox[{"y", "[", "2", "]"}], "^", "x2exp"}], " ", 
        RowBox[{
         RowBox[{"y", "[", "3", "]"}], "^", "x3exp"}]}]}], ";", " ", 
      RowBox[{"(*", " ", 
       RowBox[{
       "strip", " ", "the", " ", "coefficient", " ", "from", " ", "the", " ", 
        "term"}], " ", "*)"}], "\[IndentingNewLine]", 
      RowBox[{"coefficient", " ", "=", " ", 
       RowBox[{"term", "/", "type"}]}], ";", " ", 
      RowBox[{"(*", " ", 
       RowBox[{"only", " ", "the", " ", "coefficient", " ", "remains"}], " ", 
       "*)"}], "\[IndentingNewLine]", "\[IndentingNewLine]", 
      RowBox[{"list", " ", "=", " ", 
       RowBox[{"{", "}"}]}], ";", "\[IndentingNewLine]", 
      RowBox[{"For", "[", 
       RowBox[{
        RowBox[{"ii", "=", "1"}], ",", 
        RowBox[{"ii", "\[LessEqual]", 
         RowBox[{"Length", "[", "exponents", "]"}]}], ",", 
        RowBox[{"ii", "++"}], ",", "\[IndentingNewLine]", 
        RowBox[{
         RowBox[{"list", " ", "=", "  ", 
          RowBox[{"Join", "[", 
           RowBox[{"list", ",", 
            RowBox[{"ConstantArray", "[", 
             RowBox[{"ii", ",", 
              RowBox[{"exponents", "[", 
               RowBox[{"[", "ii", "]"}], "]"}]}], "]"}]}], "]"}]}], ";"}]}], 
       "\[IndentingNewLine]", "]"}], ";", "\[IndentingNewLine]", 
      RowBox[{"orders", " ", "=", " ", 
       RowBox[{"Permutations", "[", "list", "]"}]}], ";", 
      RowBox[{"(*", " ", 
       RowBox[{
       "Get", " ", "permutations", " ", "of", " ", "the", " ", "term"}], " ", 
       "*)"}], "\[IndentingNewLine]", 
      RowBox[{"divisor", " ", "=", " ", 
       RowBox[{"Length", "[", "orders", "]"}]}], ";", " ", 
      RowBox[{"(*", " ", 
       RowBox[{"count", " ", "number", " ", "of", " ", "permutations"}], " ", 
       "*)"}], "\[IndentingNewLine]", 
      RowBox[{"permutations", "=", "0"}], ";", "\[IndentingNewLine]", 
      RowBox[{"For", "[", 
       RowBox[{
        RowBox[{"j", "=", "1"}], ",", 
        RowBox[{"j", "\[LessEqual]", " ", 
         RowBox[{"Length", "[", "orders", "]"}]}], ",", 
        RowBox[{"j", "++"}], ",", "\[IndentingNewLine]", 
        RowBox[{
         RowBox[{"order", "=", 
          RowBox[{"orders", "[", 
           RowBox[{"[", "j", "]"}], "]"}]}], ";", "\[IndentingNewLine]", 
         RowBox[{"current", " ", "=", " ", "id"}], ";", "\[IndentingNewLine]", 
         RowBox[{"For", "[", 
          RowBox[{
           RowBox[{"ii", "=", "1"}], ",", 
           RowBox[{"ii", "\[LessEqual]", "t"}], ",", 
           RowBox[{"ii", "++"}], ",", "\[IndentingNewLine]", 
           RowBox[{
            RowBox[{"current", " ", "=", " ", 
             RowBox[{"current", ".", 
              RowBox[{"x", "[", 
               RowBox[{"order", "[", 
                RowBox[{"[", "ii", "]"}], "]"}], "]"}]}]}], ";"}]}], " ", 
          RowBox[{"(*", " ", 
           RowBox[{
           "add", " ", "each", " ", "observation", " ", "to", " ", "the", " ",
             "product", " ", "in", " ", "order"}], " ", "*)"}], 
          "\[IndentingNewLine]", "]"}], ";", "\[IndentingNewLine]", 
         RowBox[{"permutations", " ", "+=", " ", "current"}], ";"}]}], " ", 
       RowBox[{"(*", " ", 
        RowBox[{"sum", " ", "all", " ", "permutations"}], " ", "*)"}], 
       "\[IndentingNewLine]", "]"}], ";", "\[IndentingNewLine]", 
      RowBox[{"term", " ", "=", " ", 
       RowBox[{"coefficient", " ", 
        RowBox[{"permutations", " ", "/", 
         RowBox[{"(", "divisor", ")"}]}]}]}], ";", " ", 
      RowBox[{"(*", " ", 
       RowBox[{"add", " ", "coefficient", " ", "and", " ", "average"}], " ", 
       "*)"}], "\[IndentingNewLine]", 
      RowBox[{"P", "+=", "term"}], ";"}]}], "\[IndentingNewLine]", "]"}], 
   ";"}], "\[IndentingNewLine]"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"endTime", "=", 
   RowBox[{"AbsoluteTime", "[", "]"}]}], ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{
   RowBox[{"Print", "[", 
    RowBox[{"{", 
     RowBox[{"t", ",", 
      RowBox[{"endTime", "-", "startTime"}]}], "}"}], "]"}], ";"}], " ", 
  "\[IndentingNewLine]", 
  RowBox[{"(*", " ", 
   RowBox[{"implementation", " ", "of", " ", "old", " ", "formula"}], " ", 
   "*)"}]}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"startTime", "=", 
   RowBox[{"AbsoluteTime", "[", "]"}]}], ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"d0", " ", "=", " ", 
   RowBox[{
    RowBox[{"(", 
     RowBox[{
      RowBox[{"x", "[", "1", "]"}], "+", 
      RowBox[{"x", "[", "2", "]"}], "+", 
      RowBox[{"x", "[", "3", "]"}]}], ")"}], "/", "3"}]}], 
  ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"a1", " ", "=", " ", 
   RowBox[{"1", "/", 
    RowBox[{"Sqrt", "[", 
     RowBox[{"2", "k"}], "]"}]}]}], ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{
   RowBox[{"a2", " ", "=", " ", 
    RowBox[{"-", "a1"}]}], ";"}], 
  "\[IndentingNewLine]"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"xxx1", " ", "=", " ", 
   RowBox[{"MatrixPower", "[", 
    RowBox[{
     RowBox[{"(", 
      RowBox[{"d0", " ", "+", " ", 
       RowBox[{"\[ImaginaryI]", 
        RowBox[{"(", 
         RowBox[{
          RowBox[{"a1", " ", 
           RowBox[{"x", "[", "1", "]"}]}], " ", "+", " ", 
          RowBox[{"a2", " ", 
           RowBox[{"x", "[", "2", "]"}]}]}], ")"}]}]}], ")"}], ",", "t"}], 
    "]"}]}], " ", ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"xxx2", " ", "=", " ", 
   RowBox[{"MatrixPower", "[", 
    RowBox[{
     RowBox[{"(", 
      RowBox[{"d0", " ", "+", " ", 
       RowBox[{"\[ImaginaryI]", 
        RowBox[{"(", 
         RowBox[{
          RowBox[{"a1", " ", 
           RowBox[{"x", "[", "1", "]"}]}], " ", "+", " ", 
          RowBox[{"a2", " ", 
           RowBox[{"x", "[", "3", "]"}]}]}], ")"}]}]}], ")"}], ",", "t"}], 
    "]"}]}], " ", ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"xxx3", " ", "=", " ", 
   RowBox[{"MatrixPower", "[", 
    RowBox[{
     RowBox[{"(", 
      RowBox[{"d0", " ", "+", " ", 
       RowBox[{"\[ImaginaryI]", 
        RowBox[{"(", 
         RowBox[{
          RowBox[{"a1", " ", 
           RowBox[{"x", "[", "2", "]"}]}], " ", "+", " ", 
          RowBox[{"a2", " ", 
           RowBox[{"x", "[", "3", "]"}]}]}], ")"}]}]}], ")"}], ",", "t"}], 
    "]"}]}], " ", ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"P2", " ", "=", " ", 
   RowBox[{"ComplexExpand", "[", 
    RowBox[{"(", 
     RowBox[{
      RowBox[{"Re", "[", 
       RowBox[{"(", 
        RowBox[{"xxx1", "+", "xxx2", "+", "xxx3"}], ")"}], "]"}], "/", "3"}], 
     ")"}], "]"}]}], ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"endTime", "=", 
   RowBox[{"AbsoluteTime", "[", "]"}]}], ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{
   RowBox[{"Print", "[", 
    RowBox[{"{", 
     RowBox[{"t", ",", 
      RowBox[{"endTime", "-", "startTime"}]}], "}"}], "]"}], ";"}], " ", 
  "\[IndentingNewLine]", 
  RowBox[{"(*", " ", 
   RowBox[{
   "calculate", " ", "covariance", " ", "matrix", " ", "based", " ", "on", 
    " ", "new", " ", "formula"}], " ", "*)"}]}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"startTime", " ", "=", " ", 
   RowBox[{"AbsoluteTime", "[", "]"}]}], ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"vector1", "=", " ", 
   RowBox[{"Flatten", "[", 
    RowBox[{"(", 
     RowBox[{"P", "-", 
      RowBox[{"MatrixPower", "[", 
       RowBox[{"m", ",", "t"}], "]"}]}], ")"}], "]"}]}], 
  ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"varvec1", " ", "=", " ", 
   RowBox[{
    RowBox[{"Partition", "[", 
     RowBox[{"vector1", ",", "1"}], "]"}], ".", 
    RowBox[{"Partition", "[", 
     RowBox[{"vector1", ",", "4"}], "]"}]}]}], ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"way1", "=", 
   RowBox[{"Table", "[", 
    RowBox[{
     RowBox[{"Expectation", "[", 
      RowBox[{
       RowBox[{"Expectation", "[", 
        RowBox[{
         RowBox[{"Expectation", "[", 
          RowBox[{
           RowBox[{"Expand", "[", 
            RowBox[{"varvec1", "[", 
             RowBox[{"[", 
              RowBox[{"i", ",", "j"}], "]"}], "]"}], "]"}], ",", 
           RowBox[{
            RowBox[{"{", 
             RowBox[{"x111", ",", "x112", ",", "x121", ",", "x122"}], "}"}], 
            "\[Distributed]", 
            RowBox[{"MultinormalDistribution", "[", 
             RowBox[{
              RowBox[{"{", 
               RowBox[{"m1", ",", "m2", ",", "m3", ",", "m4"}], "}"}], ",", 
              RowBox[{
               SuperscriptBox["\[Sigma]", "2"], 
               RowBox[{"IdentityMatrix", "[", "4", "]"}]}]}], "]"}]}]}], 
          "]"}], ",", 
         RowBox[{
          RowBox[{"{", 
           RowBox[{"x211", ",", "x212", ",", "x221", ",", "x222"}], "}"}], 
          "\[Distributed]", 
          RowBox[{"MultinormalDistribution", "[", 
           RowBox[{
            RowBox[{"{", 
             RowBox[{"m1", ",", "m2", ",", "m3", ",", "m4"}], "}"}], ",", 
            RowBox[{
             SuperscriptBox["\[Sigma]", "2"], 
             RowBox[{"IdentityMatrix", "[", "4", "]"}]}]}], "]"}]}]}], "]"}], 
       ",", 
       RowBox[{
        RowBox[{"{", 
         RowBox[{"x311", ",", "x312", ",", "x321", ",", "x322"}], "}"}], 
        "\[Distributed]", 
        RowBox[{"MultinormalDistribution", "[", 
         RowBox[{
          RowBox[{"{", 
           RowBox[{"m1", ",", "m2", ",", "m3", ",", "m4"}], "}"}], ",", 
          RowBox[{
           SuperscriptBox["\[Sigma]", "2"], 
           RowBox[{"IdentityMatrix", "[", "4", "]"}]}]}], "]"}]}]}], "]"}], 
     ",", 
     RowBox[{"{", 
      RowBox[{"i", ",", "1", ",", "4"}], "}"}], ",", 
     RowBox[{"{", 
      RowBox[{"j", ",", "1", ",", "4"}], "}"}]}], "]"}]}], 
  ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"endTime", "=", 
   RowBox[{"AbsoluteTime", "[", "]"}]}], ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{
   RowBox[{"Print", "[", 
    RowBox[{"{", 
     RowBox[{"t", ",", 
      RowBox[{"endTime", "-", "startTime"}]}], "}"}], "]"}], ";"}], " ", 
  "\[IndentingNewLine]", 
  RowBox[{"(*", " ", 
   RowBox[{
   "calculate", " ", "covariance", " ", "matrix", " ", "based", " ", "on", 
    " ", "old", " ", "formula"}], " ", "*)"}]}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"startTime", " ", "=", " ", 
   RowBox[{"AbsoluteTime", "[", "]"}]}], ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"vector2", "=", " ", 
   RowBox[{"Flatten", "[", 
    RowBox[{"(", 
     RowBox[{"P2", "-", 
      RowBox[{"MatrixPower", "[", 
       RowBox[{"m", ",", "t"}], "]"}]}], ")"}], "]"}]}], 
  ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"varvec2", " ", "=", " ", 
   RowBox[{
    RowBox[{"Partition", "[", 
     RowBox[{"vector2", ",", "1"}], "]"}], ".", 
    RowBox[{"Partition", "[", 
     RowBox[{"vector2", ",", "4"}], "]"}]}]}], ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"way2", "=", 
   RowBox[{"Table", "[", 
    RowBox[{
     RowBox[{"Expectation", "[", 
      RowBox[{
       RowBox[{"Expectation", "[", 
        RowBox[{
         RowBox[{"Expectation", "[", 
          RowBox[{
           RowBox[{"Expand", "[", 
            RowBox[{"varvec2", "[", 
             RowBox[{"[", 
              RowBox[{"i", ",", "j"}], "]"}], "]"}], "]"}], ",", 
           RowBox[{
            RowBox[{"{", 
             RowBox[{"x111", ",", "x112", ",", "x121", ",", "x122"}], "}"}], 
            "\[Distributed]", 
            RowBox[{"MultinormalDistribution", "[", 
             RowBox[{
              RowBox[{"{", 
               RowBox[{"m1", ",", "m2", ",", "m3", ",", "m4"}], "}"}], ",", 
              RowBox[{
               SuperscriptBox["\[Sigma]", "2"], 
               RowBox[{"IdentityMatrix", "[", "4", "]"}]}]}], "]"}]}]}], 
          "]"}], ",", 
         RowBox[{
          RowBox[{"{", 
           RowBox[{"x211", ",", "x212", ",", "x221", ",", "x222"}], "}"}], 
          "\[Distributed]", 
          RowBox[{"MultinormalDistribution", "[", 
           RowBox[{
            RowBox[{"{", 
             RowBox[{"m1", ",", "m2", ",", "m3", ",", "m4"}], "}"}], ",", 
            RowBox[{
             SuperscriptBox["\[Sigma]", "2"], 
             RowBox[{"IdentityMatrix", "[", "4", "]"}]}]}], "]"}]}]}], "]"}], 
       ",", 
       RowBox[{
        RowBox[{"{", 
         RowBox[{"x311", ",", "x312", ",", "x321", ",", "x322"}], "}"}], 
        "\[Distributed]", 
        RowBox[{"MultinormalDistribution", "[", 
         RowBox[{
          RowBox[{"{", 
           RowBox[{"m1", ",", "m2", ",", "m3", ",", "m4"}], "}"}], ",", 
          RowBox[{
           SuperscriptBox["\[Sigma]", "2"], 
           RowBox[{"IdentityMatrix", "[", "4", "]"}]}]}], "]"}]}]}], "]"}], 
     ",", 
     RowBox[{"{", 
      RowBox[{"i", ",", "1", ",", "4"}], "}"}], ",", 
     RowBox[{"{", 
      RowBox[{"j", ",", "1", ",", "4"}], "}"}]}], "]"}]}], 
  ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"endTime", "=", 
   RowBox[{"AbsoluteTime", "[", "]"}]}], ";"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{
   RowBox[{"Print", "[", 
    RowBox[{"{", 
     RowBox[{"t", ",", 
      RowBox[{"endTime", "-", "startTime"}]}], "}"}], "]"}], ";"}], " ", 
  "\[IndentingNewLine]", 
  RowBox[{"(*", " ", 
   RowBox[{"compare", " ", "new", " ", "vs", " ", "old"}], " ", 
   "*)"}]}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"diff", " ", "=", " ", 
   RowBox[{"way2", "-", "way1"}]}], ";"}], "\[IndentingNewLine]", 
 RowBox[{"Eigenvalues", "[", "diff", "]"}], "\[IndentingNewLine]", 
 RowBox[{
  RowBox[{"obj", "=", 
   RowBox[{"EvaluationCell", "[", "]"}]}], ";"}]}], "Input",
 CellChangeTimes->{
  3.694366166519425*^9, {3.694366731629274*^9, 3.694366774487414*^9}, {
   3.6943672962346725`*^9, 3.6943672970608463`*^9}, {3.6943686944237223`*^9, 
   3.6943686985162296`*^9}, {3.694394447390716*^9, 3.694394447484427*^9}, 
   3.6974053752941437`*^9, 3.697497898767414*^9, 3.6974980487382574`*^9, {
   3.697500999985654*^9, 3.6975010524740953`*^9}, {3.697501089778041*^9, 
   3.697501132097255*^9}, {3.697501477703276*^9, 3.6975015190992317`*^9}}],

Cell[CellGroupData[{

Cell[BoxData[
 RowBox[{"{", 
  RowBox[{"2", ",", "0.0469107`6.122816907020173"}], "}"}]], "Print",
 CellChangeTimes->{3.6975015215244274`*^9}],

Cell[BoxData[
 RowBox[{"{", 
  RowBox[{"2", ",", "0.0312516`5.946417250484317"}], "}"}]], "Print",
 CellChangeTimes->{3.6975015215556707`*^9}],

Cell[BoxData[
 RowBox[{"{", 
  RowBox[{"2", ",", "0.3561947`7.003232446585353"}], "}"}]], "Print",
 CellChangeTimes->{3.6975015218962755`*^9}],

Cell[BoxData[
 RowBox[{"{", 
  RowBox[{"2", ",", "0.2968801`6.924126081122877"}], "}"}]], "Print",
 CellChangeTimes->{3.6975015221931553`*^9}]
}, Open  ]]
}, Open  ]],

Cell[CellGroupData[{

Cell[BoxData[
 RowBox[{"CellToTeX", "[", 
  RowBox[{
   RowBox[{"NotebookRead", "[", "obj", "]"}], ",", 
   RowBox[{"\"\<Style\>\"", "\[Rule]", "\"\<Input\>\""}]}], "]"}]], "Input",
 CellChangeTimes->{{3.69749797902544*^9, 3.697497980270496*^9}, {
  3.6975015794497375`*^9, 3.69750158073456*^9}}],

Cell[BoxData["\<\"\\\\begin{mmaCell}[moredefined={k, t, s1, i, s2, x, x111, \
x112, x121, x122, x211, x212, x221, x222, x311, x312, x321, x322, m, id, \
startTime, P, terms, j, term, x1exp, x2exp, x3exp, exponents, type, \
coefficient, list, ii, orders, divisor, permutations, order, current, \
endTime, d0, a1, a2, xxx1, xxx2, xxx3, P2, vector1, varvec1, way1, vector2, \
varvec2, way2, diff, obj}]{Input}\\n  k=3;\\n  t = 2;\\n  s1 = \
\\\\mmaUnderOver{\\\\(\\\\pmb{\\\\sum}\\\\)}{i=1}{k}(y[i])/k;\\n  s2 = \
\\\\mmaUnderOver{\\\\(\\\\pmb{\\\\sum}\\\\)}{i=1}{k}(y[i]-s1)^2/k;\\n  (* 2x2 \
matrix *)\\n  x[1] = (x111\\tx112\\n  x121\\tx122);\\n  x[2] = \
(x211\\tx212\\n  x221\\tx222);\\n  x[3] = (x311\\tx312\\n  x321\\tx322);\\n  \
m=(m1\\tm2\\n  m3\\tm4);\\n  id = (1\\t0\\n  0\\t1);\\n  (* implementation of \
new formula *)\\n  startTime=AbsoluteTime[];\\n  P=0;\\n  terms = \
MonomialList[Sum[\\\\mmaFrac{\\\\mmaSup{(-1)}{i} \
(2i-1)!!Binomial[t,2i]}{Product[k+2j-3,\\\\{j,1,i\\\\}]}\\\\mmaSup{s1}{t-2i}\\\
\\mmaSup{s2}{i},\\\\{i,0,Floor[t/2]\\\\}]];  (* list of monomial terms in P \
*)\\n  \\n  For[i=1,i\\\\(\\\\pmb{\\\\leq}\\\\)Length[terms],i++,\\n  term = \
terms[[i]];\\n  x1exp = Exponent[term,y[1]];\\n  x2exp = Exponent[term,y[2]];\
\\n  x3exp = Exponent[term,y[3]];\\n  \\n  exponents=\\\\{x1exp,x2exp,x3exp\\\
\\}; (* list of the exponents appearing in this type of term *)\\n  \\n  type \
= y[1]^x1exp y[2]^x2exp y[3]^x3exp; (* strip the coefficient from the term *)\
\\n  coefficient = term/type; (* only the coefficient remains *)\\n  \\n  \
list = \\\\{\\\\};\\n  \
For[ii=1,ii\\\\(\\\\pmb{\\\\leq}\\\\)Length[exponents],ii++,\\n  list =  \
Join[list,ConstantArray[ii,exponents[[ii]]]];\\n  ];\\n  orders = \
Permutations[list];(* Get permutations of the term *)\\n  divisor = \
Length[orders]; (* count number of permutations *)\\n  permutations=0;\\n  \
For[j=1,j\\\\(\\\\pmb{\\\\leq}\\\\) Length[orders],j++,\\n  \
order=orders[[j]];\\n  current = id;\\n  \
For[ii=1,ii\\\\(\\\\pmb{\\\\leq}\\\\)t,ii++,\\n  current = \
current.x[order[[ii]]]; (* add each observation to the product in order *)\\n \
 ];\\n  permutations += current; (* sum all permutations *)\\n  ];\\n  term = \
coefficient permutations /(divisor); (* add coefficient and average *)\\n  \
P+=term;\\n  ];\\n  \\n  endTime=AbsoluteTime[];\\n  \
Print[\\\\{t,endTime-startTime\\\\}]; \\n  (* implementation of old formula \
*)\\n  startTime=AbsoluteTime[];\\n  d0 = (x[1]+x[2]+x[3])/3;\\n  a1 = \
1/Sqrt[2k];\\n  a2 = -a1;\\n  \\n  xxx1 = MatrixPower[(d0 + \\\\mmaDef{i}(a1 \
x[1] + a2 x[2])),t] ;\\n  xxx2 = MatrixPower[(d0 + \\\\mmaDef{i}(a1 x[1] + a2 \
x[3])),t] ;\\n  xxx3 = MatrixPower[(d0 + \\\\mmaDef{i}(a1 x[2] + a2 x[3])),t] \
;\\n  P2 = ComplexExpand[(Re[(xxx1+xxx2+xxx3)]/3)];\\n  \
endTime=AbsoluteTime[];\\n  Print[\\\\{t,endTime-startTime\\\\}]; \\n  (* \
calculate covariance matrix based on new formula *)\\n  startTime = \
AbsoluteTime[];\\n  vector1= Flatten[(P-MatrixPower[m,t])];\\n  varvec1 = \
Partition[vector1,1].Partition[vector1,4];\\n  \
way1=Table[Expectation[Expectation[Expectation[Expand[varvec1[[i,j]]],\\\\{\
x111,x112,x121,x122\\\\}\\\\(\\\\pmb{\\\\unicode{f3d2}}\\\\)\
MultinormalDistribution[\\\\{m1,m2,m3,m4\\\\},\\\\mmaSup{\\\\mmaUnd{\\\\(\\\\\
pmb{\\\\sigma}\\\\)}}{2}IdentityMatrix[4]]],\\\\{x211,x212,x221,x222\\\\}\\\\(\
\\\\pmb{\\\\unicode{f3d2}}\\\\)MultinormalDistribution[\\\\{m1,m2,m3,m4\\\\},\
\\\\mmaSup{\\\\mmaUnd{\\\\(\\\\pmb{\\\\sigma}\\\\)}}{2}IdentityMatrix[4]]],\\\
\\{x311,x312,x321,x322\\\\}\\\\(\\\\pmb{\\\\unicode{f3d2}}\\\\)\
MultinormalDistribution[\\\\{m1,m2,m3,m4\\\\},\\\\mmaSup{\\\\mmaUnd{\\\\(\\\\\
pmb{\\\\sigma}\\\\)}}{2}IdentityMatrix[4]]],\\\\{i,1,4\\\\},\\\\{j,1,4\\\\}];\
\\n  endTime=AbsoluteTime[];\\n  Print[\\\\{t,endTime-startTime\\\\}]; \\n  \
(* calculate covariance matrix based on old formula *)\\n  startTime = \
AbsoluteTime[];\\n  vector2= Flatten[(P2-MatrixPower[m,t])];\\n  varvec2 = \
Partition[vector2,1].Partition[vector2,4];\\n  \
way2=Table[Expectation[Expectation[Expectation[Expand[varvec2[[i,j]]],\\\\{\
x111,x112,x121,x122\\\\}\\\\(\\\\pmb{\\\\unicode{f3d2}}\\\\)\
MultinormalDistribution[\\\\{m1,m2,m3,m4\\\\},\\\\mmaSup{\\\\mmaUnd{\\\\(\\\\\
pmb{\\\\sigma}\\\\)}}{2}IdentityMatrix[4]]],\\\\{x211,x212,x221,x222\\\\}\\\\(\
\\\\pmb{\\\\unicode{f3d2}}\\\\)MultinormalDistribution[\\\\{m1,m2,m3,m4\\\\},\
\\\\mmaSup{\\\\mmaUnd{\\\\(\\\\pmb{\\\\sigma}\\\\)}}{2}IdentityMatrix[4]]],\\\
\\{x311,x312,x321,x322\\\\}\\\\(\\\\pmb{\\\\unicode{f3d2}}\\\\)\
MultinormalDistribution[\\\\{m1,m2,m3,m4\\\\},\\\\mmaSup{\\\\mmaUnd{\\\\(\\\\\
pmb{\\\\sigma}\\\\)}}{2}IdentityMatrix[4]]],\\\\{i,1,4\\\\},\\\\{j,1,4\\\\}];\
\\n  endTime=AbsoluteTime[];\\n  Print[\\\\{t,endTime-startTime\\\\}]; \\n  \
(* compare new vs old *)\\n  diff = way2-way1;\\n  Eigenvalues[diff]\\n  \
obj=EvaluationCell[];\\n\\\\end{mmaCell}\"\>"], "Output",
 CellChangeTimes->{3.6975015850821247`*^9}]
}, Open  ]]
},
WindowSize->{979, 666},
WindowMargins->{{142, Automatic}, {-154, Automatic}},
FrontEndVersion->"10.0 for Microsoft Windows (64-bit) (December 4, 2014)",
StyleDefinitions->"Default.nb"
]
(* End of Notebook Content *)

(* Internal cache information *)
(*CellTagsOutline
CellTagsIndex->{}
*)
(*CellTagsIndex
CellTagsIndex->{}
*)
(*NotebookFileOutline
Notebook[{
Cell[558, 20, 91, 2, 31, "Input"],
Cell[CellGroupData[{
Cell[674, 26, 18687, 511, 1982, "Input"],
Cell[CellGroupData[{
Cell[19386, 541, 142, 3, 23, "Print"],
Cell[19531, 546, 142, 3, 23, "Print"],
Cell[19676, 551, 142, 3, 23, "Print"],
Cell[19821, 556, 142, 3, 23, "Print"]
}, Open  ]]
}, Open  ]],
Cell[CellGroupData[{
Cell[20012, 565, 296, 6, 31, "Input"],
Cell[20311, 573, 4945, 68, 1992, "Output"]
}, Open  ]]
}
]
*)

(* End of internal cache information *)